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SECTION 1 


INTRODUCTION 


1.1 Purpose of Study 

In recent years, with accelerated developments in aerospace technology, 
the needs for accurate and rapid methods for flight-vehicle-safety analysis 
and for structural design have been keenly felt. In design, one must account 
for the interaction of the structural system witn various types of environ- 
mental dynamic conditions such as landing, maneuvering, gust, blast, etc. Also 
of interest is the interaction of a high-velocity fragment with the structure 
which is intended to contain and/or to deflect the fragment; uncontained frag- 
ments wliicn may be generated from the failure of high-speed rotating turbojet 
aircraft engine parts could damage equipment and threaten passenger safety [1]*. 

For efficient minimum weight (optimum) design, it is often necessary to 
take full advantage of the load-carrying capacity of available materials by 
permitting the material to proceed well into the plastic range; thus, non- 
linear material behavior must be accounted for in some design analyses. Also, 
frequently, it is necessary in predicting structural responses to consider geo- 
metric nonlinearities, specifically, when the deflections are large enough. 
Multilayer structures involving composite materials which have high strength and 
low density properties are becoming useful in aerospace structures as well as in 
other applications. The complex and nonlinear character of such structural 
problems, however, makes it almost impossible for one to carry out a conventional 
analytical solution in closed form. In order to provide a practical means of 
obtaining meaningful predictions for this kind of complex problem, numerical 
analysis procedures have been developed. 

The general numerical methods of structural analysis may be divided con- 
veniently into two categories. The first category, "numerical solution of the 
governing algebraic and/or differential equations" is based on mathematically 


★ 

Numbers in brackets [ ] refer to references cited at the end of the text. 
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approximating and solving the governing differential equations either by 
finite-differences [2, 3, for example] or by forward integration* 14, 5, 6, for 
example] . The second category is the "finite-element method" which is most 
systematically based on variational principles. In this method, the solid con- 
tinuum is represented by a finite number of regions which are connected at a 
finite number of nodes along interelement boundaries ; the geometric and the 
material properties of the continuum may be faithfully retained in the idealized 
structural assembly. In the past several years , the finite-element method of 
analysis has undergone intensive development and has proved to be a very effec- 
tive and powerful method for analyzing certain classes of problems, especially 
for a continuum or structure with complicated boundary conditions, geometric 
shape, and material properties. The relative ease and versatility with which 
the finite element (FE) method can be applied to complex structural shapes in 
comparison with the finite-difference (FD) method is often regarded as an im- 
portant attribute of the finite-element method. Accordingly, the assumed-dis- 
placement version of the finite-element method of analysis for large-deflection 
elastic-plastic transient responses of simple structures is developed in the 
present study*. Various aspects of the present development are given in 
Sections 2 through 5; for convenience, a review of the pertinent literature on 
this topic is given in Subsection 3.1. 

Discussed in Section 6 is a timely problem to which the present method 
of analysis has been applied; this concerns rotor-blade fragment impingement 
upon (a) a complete circular structural containment ring or (b) a segment of a 
circular structural ring which is supported in one of several ways. Sought 
for (a) is a prediction of the motion and structural deformation of the contain- 
ment ring, as well as the motion of the (idealized-as-rigid in the present 
first-approximation analysis) rotor blade; suitably accurate engineering pre- 
dictions of strains are desired. Similar data in category (b) are sought, but 
there is greater interest in observing the "diverted trajectory" of the rotor 

+ This method is applicable to problems in one space variable such as beams and 
shells of revolution. 

* 

Principal emphasis has been devoted to employing the assumed-displacement 
type of finite-element model; however, other types of finite-element models 
such as equilibrium-stress models, hybrid models, and mixed models [7] 
could also be used for this type of problem. 
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blade — the idea of interest being to divert the blade to a new "harmless 
direction" . Further background on the rotor fragment containment/deflection 
problem and on alternate or supplementary methods for its analysis is given 
in Subsection 6.1. Section 6 is devoted to discussing the rotor burst frag- 
ment containment/deflection problem, and the application of a method of analy- 
sis developed in this study to analyze some simple cases of fragment impinge- 
ment upon containment and/or deflection structures, 

1.2 Synopsis of Investigation 

The present study is devoted to developing and validating a method to 
analyze the large deflection transient responses of simple structures , in- 
cluding elastic-plastic material behavior . The assumed-displacement finite- 
element approach which is based upon the Principle of Virtual Work and 
D'Alembert's Principle (or equivalently, Hamilton's Principle) is used to formu- 
late the governing equations. The resulting equations of motion are solved by 
a direct timewise numerical integration scheme. This method has been evaluated 
(see Section 5) by applying it to a sequence of transient response examples 
having reliable independent analytical or experimental results, for a definitive 
comparison. This development is then applied to the rotor fragment containment/ 
deflection problem. A simplified fragment structure interaction model which 
utilizes momentum and energy considerations of the system is employed to pre- 
dict the collision-induced velocities imparted "instantaneously” to the 
affected ring segment and to the fragment. 

Sections 2 through 5 pertain to the development and evaluation of the 
present finite-element method for predicting both small-deflection and large- 
deflection elastic and/or elastic-plastic transient responses of simple struc- 
tures. The general equations which govern large elastic-plastic dynamically- 
induced deformations of a continuum are presented in Section 2 in general three- 
dimensional tensor form. Section 3 is devoted to the development of an over- 
all method of solution, following a review of the pertinent literature dealing 
with this class of problems; the spatial finite-element approximation together 
with the temporal finite-difference approximation are used. Also in Section 3, 
the equations of motion for the finite-element treatment are derived from a 
variational statement consisting of the Principle of Virtual Work and 
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D'Alembert's Principle; the resulting equations are developed in two forms: 

(a) the "conventional form" (see Eq. 3.1) , and (b) an "improved" unconventional 
form (see Eq. 3.2). The new improved formulation (b) is shown (in Section 5) 
to be more efficient and economical for computing the large-deflection elastic- 
plastic transient responses of simple structures than is the conventional 
finite-element formulation (a). In Section 4, the general formulation dis- 
cussed in Section 3 is developed in detail for an arbitrarily curved beam 
having (a) zero and (b) non-zero transverse shear deformation. Finally, 

Section 5 contains an assessment of this method of analysis by means of a 
sequence of problems for beam and ring example structures which are subjected 
to transient mechanical loading or to initial impulsive loading; comparisons 
of the present predictions are made with reliable experimental data and/or 
independent predictions (finite-difference and/or analytical) . 

In Section 6, the problem of burst- rotor fragments interacting with 
either a complete or a partial containment/deflection structural ring is dis- 
cussed. Energy and momentum considerations are employed to predict the collision- 
induced velocities which are imparted to the colliding fragment and to the af- 
fected ring segment (the associated analysis method is termed the collision- 
imparted velocity method, CIVM) . This collision analysis is combined with the 
earlier-developed FE analysis to permit one to predict the resulting large de- 
formation responses of containment/deflection rings. Comparisons with limited 
experimental data are also given. 

The entire study is summarized and pertinent conclusions are drawn in 
Section 7. 

Also, three appendices are included. Appendix A contains the description 
of the mechanical sublayer model for strain-hardening, strain-rate sensitive 
material behavior. In Appendix B, all of the matrices used in the present 
analysis for the assumed displacement finite-element formulations are listed. 
Appendix C contains, for purposes of illustration, the formulation of the 
equations of motion based on a mixed finite-element method (1) by assuming a 
displacement field which is continuous over the entire solid and (2) by using 
an assumed-stress field for the individual element; however, no evaluation of 
this model has been made in this report. 
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SECTION 2 


GENERAL FORMULATION 


2.1 Introduction 

In this section, the equations which govern the large deflection dynamic 
and/or static responses of a continuum are presented. Elastic-plastic strain- 
rate-dependent material behavior is considered, but not viscoelasticity. The 
equations are derived in general tensor form for convenience and generality 
and, therefore, any coordinate system can be employed for defining the three- 
dimensional space containing the continuum. These equations are later speci- 
alized to treat simpler classes of problems. 

The terminology "large deflections" sis used here indicates, for example, 
that the lateral deflections of beams and/or plates are large compared with 
the thickness of the structure; the change in geometry is significant. Through 
the strain-displacement relations and the equilibrium equations, the geometric 
nonlinearities are introduced into the theory. It should be noted, however, 
that in the present analysis, the strains (extension and shear) are treated as 
being small compared with unity. 

In a finite-strain analysis, there are various possible types of defini- 
tions of stresses and strains (Refs. 8, 9, 10) based on either the predeformed 
or post-deformed configurations of the continuum, and the distinctions between 
them cannot be neglected; whereas, they may be indistinguishable in small- 
strain theory. The constitutive relations for cases involving finite elastic- 
plastic strain are uncertain (not well established) , primarily because of the 
lack of adequate experimental data. However, finite-strain predictions have 
been of increasing interest recently, for instance, in connection with ex- 
plosive forming of structural shapes. The analysis tools for this class of 
problem are much desired, and further research is required. 

In this study, the indicial notation and summation convention associated 
with vector and tensor analysis are used. 
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2. 2 Governing Equations 


2.2.1 The Strain-Displacement Relations 


Consider the continuum in its process of deformation from an initially 
undeformed state which exists at time t (Fig. '1). Each material point in the 
3-dimensional space can be identified by a general curvilinear coordinate sys- 
tem (termed Lagrangian) and an inertial Cartesian rectangular coordinate 
system fixed in space. Let the position vector of a given material point 

in the initial undeformed state be r (£ 3 , t ) and the position vector of the 

o 

same material point in the deformed* state at any instant of time t be 
R(S j , t). 

Base vectors for the coordinate system may be defined in both the un- 
deformed state and the deformed state of the continuum, respectively, as 


h- 


9 r 
d ** 


and 


pr = 13 


( 2 . 1 ) 


The metric tensor g. . associated with the undeformed state and the metric 
tensor G ±j associated with the deformed state are defined, respectively, by: 


fij ~ fx ff 


% " 5* • 


( 2 . 2 ) 


mn inn 

The contravariant metric tensors g and G are defined, respectively, through 


the relations** 

mn r h 

(tmf ft 


tnn r n 

= h 


(2.3) 


where 6 is the Kronecker delta defined as 
P 


^ = 


if n * p 
if n + p 


(2.4) 


The respective contravariant base vectors are defined by 

j* = f’¥ r , G = J 5, (2.5) 

* 

Quantities associated with the deformed state are represented by capital 
letters, and lower case letters denote the initial undeformed state. 

** 

Index quantities such as i, j, p, m, n, k, etc., take on values 1, 2, and 3. 
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The deformed position vector R and the undeformed position vector r 


are related by the displacement vector t) as follows: 

r (^, t) = 0 ) + t) 


( 2 . 6 ) 


The displacement vector v may be expressed in component form with re- 


spect to the base vector system of the undeformed state at time t as 


v = * n 


- Vj 


(2.7) 


The Lagrangian strain tensor associated with the deformation process 
is defined by 

~Yxf ~ T ( fct ) 


(2.3) 


and can be written in terms of the displacement components as 

^ V p* + V ’*■ V k>j ) (2.9) 

where v. , denotes covariant differentiation of v. with respect to Ep using the 
metric tensor of the undeformed state and is defined as 



( 2 . 10 ) 


where {i j} is the Christoffel symbol of the second kind associated with the 
o 

undeformed state (subscript "o") and is defined by 



( 2 . 11 ) 


Differentiating Eqs. 2.7 and 2.9 with respect to time t which is a monotonically 
increasing quantity, and using the fact that g^ and {i j are not functions of 
time (fixed quantities with respect to the initial state) , one may obtain the 
following expression for the displacement rate: 


V 




( 2 . 12 ) 


Hence, the strain rate y. . which is related to the displacement and displaee- 

13 

ment rate is given by 
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-v - -L f ( (S' 1 *5 + V k.) V . +( fik 4 v^ ) v, . 1 

rxf ~ i O k,f lu j >y k,+ ) 


1 


(2.13) 


where ( } denotes partial differentiation with respect to time t. 


In the problem of interest, expressions are required for the strain incre- 
ment Ay^ and the displacement increment Av which occur when the continuum is 
deformed incrementally from one deformed state at time t^ to another deformed 
state at time where = t^ + At (the solution procedure is to be dis- 
cussed later for transient-response behavior of a continuum) . If the time- 

2 — 

step At between these two instants is small, the displacement (denoted by v) 
at time t can be constructed by employing the Taylor series expansion involv- 
ing the finite increment At in time* : 


V - v 

or in component form 

2.. ,, + ) + 1 V. + 0 (*tr 


■ = V, - Vi (*t) * 


(2.14) 


(2.15) 


and correspondingly, the Lagrangian strain y^> at time t 2 may be expressed 


as 




16) 


In Eq. 2.16 terms of the or^er of (At) J and higher are neglected because of 
the smallness of the time-step At. 

It perhaps should be mentioned that the displacement vector v can be ex- 
pressed alternatively in terms of components with respect to the deformed base 
vectors at time t as 

v = V = V m Q (2.i 7 ) 

and the corresponding strain tensor (termed Almansi, see Ref. 11) has the form 


* 

Note that prescripts, both sub and super, are used here to identify quantities 
associated .with particular instants in time. 
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(2.18) 


■y„„ = J ( v m i„ + V„|„ - v fe L v k u ) 

where V . indicates also covariant differentiation but with respect to the 

m jn r 

metric tensor of the deformed state and is defined as 


V i = - 

v m| n ^ '£ j n 

J fc \ 
[mn J 

V k 

(2.19) 

f k \ d • 

U n j 2 

II 

D * 

(2.20) 


— JC 

Since G as well as { } are functions of time, the displacement vector 

m mn 

rate and the strain rate expressed in terms of components with respect to the 
deformed base vector system will not take as simple a form as Eqs. 2.12 and 
2.13 in which the components are expressed with respect to the undeformed base 
vector system. Hence, the strain-displacement relations represented by Eq. 2.9 
are chosen for use. 

2.2.2 Descriptions of Stresses and Forces on a 
Deformed Continuum Element 

Consider the geometry and the forces acting upon an arbitrary parallele- 
piped element on its path of deformation (see Figs. 1 and 2). In the initial 

undeformed state, it has the volume dV bounded by surfaces dA (having edges 
n- p- o om 

dt g n and dC P g p , where m ? n ^ p) . After deformation, this Lagrangian paral- 
lelepiped has a volume dV and is bounded by surfaces dA (having edges 
n— d— m 

d£ G , dl^G ) . Thus , 

A VC , AV = jG' 

A _ d A^jGQ^d^d^ 

( 2 . 21 ) 

where g and G are the determinants of g and G„ , respectively. 

In the following, convenient descriptions for stresses and forces acting 
upon the deformed parallelepiped element are to be discussed as a prelude for 
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their use in the Principle of Virtual Work from which the equations of motion 
and boundary conditions will be obtained. 


Let o denote the stress vector associated with the £ x surface (at 
i ° 

£ ■ const) in the deformed state measured per unit initial undeformed area . 

The Kirchhoff stress tensor S 13 is defined by (Fig. 2) : 


tJf 


-- s ; « ' 


( 2 . 22 ) 


f/JT r =s‘’(sh v *t > h 


h \ n 


(2.22a) 


— i 


The stress resultant T exerted on the deformed surfaces by the neighboring 
element of the continuum may be expressed in the form 




.t 


i, 


k + V k < ) " 


(2.23) 


Alternatively, let O be the stress vector associated with the £ surface 
in the deformed body but measured per unit deformed area; thus, the Eulerian 
stress tensor, x 13 , is defined by (Fig. 2): 

— i J ft “ — i 


Jg* = T* 7 2T, = T, s 


O' JG* 


1 -f 


Then the stress resultant, T , may be expressed as 


T* = TJGG* = JW t l> G 


(2.24) 


(2.25) 


Likewise, the body force, F (inertia, gravitational, magnetic, etc.) per unit 
mass acting on the deformed body may be described in component form either with 
respect to the undeformed base vector of the coordinate system £ X as 

F = f A % = h f ‘ (J.36) 

or with respect to the deformed base vectors of the coordinate system £ x as 

f - F X G*. = Fx G X (2.27) 

It should be noted that the stress tensors S 13 and x 13 are related to 
each other by 
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(2.28) 


G' 

t ° 

ij ij 

Also, from moment equilibrium it can be shown that S and T are symmetric 
stress tensors . Thus 




jii i *■ 

T - T 


(2.29) 


2.2.3 Conservation of Mass 


The conservation of mass may be expressed as (Ref. 11) to be 


J5T P. d ^ = ff(p 6 v 


v. v 

or using Eq. 2.21, one may write 


Wi 


?.is **:**'*? 

V 


(2.30) 


(2.31) 


where the integrals extend over the same particles. The quantity P q is the 
density in the original configuration and p is the current density of the de- 
formed body. Since this relation must hold for all bodies, one has 


Po If - P 


(2.32) 


This gives the relation of densities in different configurations of the con- 
tinuum to that of the original undeformed continuum. 


2.2.4 The Principle of Virtual Work 

In this subsection, the Principle of Virtual Work for the continuum 
will be presented and is employed to derive the equations of equilibrium and 
the boundary conditions. 

Consider the deformed continuum in equilibrium, under the action of 
body forces, externally applied surface tractions, and with arbitrary deforma- 
tion conditions consistent with the prescribed geometric boundary conditions. 
Let this equilibrium configuration be subjected to an arbitrary and independent 
set of infinitesimal virtual displacements without violating the prescribed 
geometric boundary conditions. The Principle of Virtual Work states that the 
virtual work done by the external forces (body forces and surface tractions) 
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is equal to the virtual work of the internal stresses. By utilizing the unde- 
formed metric of the coordinate system the Principle of Virtual Work may 
he stated mathematically as 


& U - S W = o 


(2.33) 


where 


SU 


/// 

V 0 


= variation of the work of 
the internal forces 


(2.33a) 


^ y\f -ff(o To' J V dA 0 ‘ variation of the work of: (1) 

i))‘° 0 J ' the body force F per unit mass 

y o Ao(r and (2) of the externally applied 

surface traction T , measured per 
o 

unit area of the undeformed state 

(2.33b) 

where S 13 is the Kirchhoff stress tensor 
■y\ j is the Lagrangian strain tensor 

6 is the usual symbol denoting the variation of a function 
6v is the variation of the displacement vector and may be 
expressed in terms of components with respect to the 
undeformed base vectors as 


$V = $(V fix) ( ^ ^ J fix (because g^ does not de- 

pend upon displacement, 

6^ - 0) (2.34) 


In Eqs. 2.33a and 2.33b, V is the entire undeformed volume of the continuum 

o 

which is bounded by the undeformed surface A q . The boundary surface A q can 
be dividad into two parts: (a) prescribed surface traction boundary A qo and 
(b) prescribed displacement boundary A Qy . Then the variation in the 
Lagrangian strain tensor, 6y. , associated with the displacement variation 6v 
about the deformed equilibrium configuration may be expressed as 
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(2.35) 



where ( ) signifies covariant differentiation with respect to the undeformed 

,n i 

metric of the coordinate system £ . 

Substituting Eq. 2.35 into Eq. 2.33a one has 

k , i/ k 


su =-{(( s‘*{ + v ,x) ^ 


(2.36) 


V. 


where the symmetry of S J has been used to combine terms. Upon integration by 
parts and using Green's theorem to convert volume integrals to surface integrals, 
Eq . 2.36 becomes 

k 

SU = 

Ao(r 

k , ,, k 


S i ^(5 t t n.f d A, 

-/)Tf 5 XJ ( s i )] .Sv k dV. 


V. 


(2.37) 


where n Q . is the component of the unit normal vector, n^, to the undeformed 
boundary surface A q , which is defined as 


n. = 

Substituting Eqs. 2.26 and 2.34 into Eq. 2.33b results in 

S W-fff pJ k Sv k dV. -rff T, k Sv k dA. 

\ f /4o<r 

v° _ 

where T q is the component of T q in the direction of the undeformed base vector. 


(2.38) 


(2.39) 
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Thus , 


T. = Tt f k = T„ k ^ 


( 2 . 40 ) 


Substituting Eqs. 2.37 and 2.39 into the Principle of Virtual Work one obtains 

+pj l ) Sv k 6V 0 

\f( Sv kd A 0 =0 


(2.41) 

Since the variation ov^ is independent and arbitrary, this leads to the follow- 
ing equation of equilibrium at any point in the continuum: 

k \ 7 . rs r k 


{S tf (Si + V ,X ] li * 


0 


(2.42) 


x i i 

which is a bilinear function in S and v . 

Also, on the boundary A^ where the surface traction is prescribed, the 
following boundary condition relation must be satisfied: 


= TJ 


k 


(2.43) 


On the portion of the boundary A qv where the displacement is prescribed, the 


relation 


x i 

V = V 


(2.44) 


should be satisfied. 


It should be noted that all pertinent quantities appearing in 6u and <$W 
of Eq. 2.33 are defined consistently with respect to the (initial) undeformed 
metric of the coordinate system S* . Also, the integrations are extended over 
the original (known) configuration of the continuum. This is usually called a 
"Lagrangian description" . 

Alternatively, by using the relations 

p.Jf = pJW, S^' = ]f fdA-fdA «.«, 
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the Principle o£ Virtual Work, Eg. 2.33, can be converted into the following 
equivalent form with reference to the present deformed state (designated as an 
Eulerian* description) : 

fffr dV=fffpF-SvdV^ff j-SvdA , 2 . 46 , 

V V A<r 

It is worthy of special mention that this basic variational formulation, 
the Principle of Virtual Work, holds independent of the material stress-strain 
relation and the existence of potential functions of the external forces. Also, 
it embodies the equation of equilibrium of the continuum (overall and/or at any 
generic point) . The variational technique has proved to be a very powerful and 
easily applied method for analyzing solid continuum problems. In the present 
study, the Principle of Virtual Work, together with the concept of D'Alembert's 
Principle is employed for analyzing the large-deflection elastic-plastic 
transient responses of structures which may be subjected to an initial impulse 
followed by a time-dependent externally-applied forcing function and geometric 
constraints . 


The D'Alembert Principle states that the dynamic system can be considered 
to be in equilibrium under the externally-applied forces if the inertial forces 
are taken into account. The Principle of Virtual Work, Eq. 2.33, now with the 
virtual work done by the inertial forces specified explicitly may be written 
as 


su -Z w I = ° 


(2.47) 


where 


ll = ((( P.R-S 7 dV. =(f(p,v-Sv AM. 

V. v. 

= variation of the work done by the inertia forces 


(2.47a) 


* 

It should be noted that, in fluid mechanics, the Eulerian description is con- 
cerned with the state of motion through a fixed grid or a fixed volume in space, 
in contrast with the above solid-mechanics terminology. 
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In Eq. 2.47a, R * v since R = r + v and r is the position vector of the 
unde formed body. It is assumed that all body forces other than inertia forces 
are taken into account within 6w. 


Let Eq. 2.47 be integrated with respect to time between two time limits 
t 1 and t 2 to represent the conventional statement of Hamilton's Principle. Thus 


{**{8 U-SW+SIJ it = o 


Integrating the last term which appears under the integral sign of Eq. 2.48 by 



(2.49) 


Since if the conditions of the dynamic system are prescribed (or are otherwise 
known) at the two time limits t^ and t 2 and are not subjected to variation, 
one has 


l 7 ( t = t, ) = 0 


S V(t=t £ )= 0 


(2.50) 



where 



= the kinetic energy of the body (as 
usually seen in Hamilton's princi- 
ple) . 


Substituting Eq. 2.51 into Eq. 2.48, one has the Principle of Virtual 
Work of the dynamic system: 
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(2.52) 


f t ‘{8U-$W- S K ] Jt = o 

It perhaps should be pointed out that either the Lagrangian form or the 
Euler ian form of Eg. 2.52 could be used in the exact or in an approximate 
analysis of structural responses. For convenient reference, these two forms 
are written out in full as follows: 

f ( /// s i y, t d -[If pj‘ 5 i v. -[( Tc s . 

1 V 0 V 0 A 0 <r 

- fff ft ** Sv x dV.] it = 0 

■ V. 

'(Lagrangian form) (2.52a) 

f‘ ( iff T ‘* s Y 'f d V ~li ? F * S Vi T ‘ S dA 

J y y Ar 

-j[J p 7-8 7 dV] dt = 0 

V 

(Eulerian form) (2.52b) 

It is seen that the Eulerian description is based on the current deformed 
volume, area, and metric tensor; these quantities are functions of the dis- 
placement vector, v, which is to be determined as a function of time. On the 
other hand, the Lagrangian description is based on the original (known) volume, 
area, and metric tensors and, accordingly , is more convenient. 

Because of the complexity in certain structural problems, numerical 
methods are of interest and include, for example, finite difference, finite 
element, forward integration, etc. In this study, however, the finite-element 
approach is used and is based on the variational formulation given by Eq. 2.52a. 
Since, in this Lagrangian point of view, the strain would take on its simple 
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form and the mass matrix (as discussed later) would be constant throughout the 
calculation of the transient response (and hence need not be re-evaluated at 
each stage in the incremental time analysis as would be required by the Eulerian 
description) , the use of the bagrangian description would appear to result in 
some computational advantages and time saving in solving the problem over the 
use of the Eulerian description. 

It perhaps should be noted that the Principle of Virtual Work, Eq. 2.52, 
can be applied either to conservative or to non-conservative dynamic systems. 

For conservative dynamic systems (that is, if the potential of the internal 
stress and externally-applied force exist) , as well as for dynamic systems with 
both conservative and non-conservative forces, another very useful energy 
principle called Hamilton's Principle could also be used; this principle can 
be derived from the Principle of Virtual Work as employed here (i.e., Eq. 2.48 
or Eq. 2.52). 

Finally, it should be noted that one need not use the time integrated 
form of the Principle of Virtual Work (within which is imbedded D'Alembert's 
Principle) if one wishes simply to obtain the correct equations of motion and 
the correct boundary conditions — the use of Eqs. 2.47 and 2.47a which hold 
at each instant in time will suffice for this purpose. However, if one 
wishes to formulate a finite-element analysis in both space and time (see 
Ref. 12), the “time integrated" variational statement represented by Eq. 2.48 
or Eq. 2.52a, for example , is very useful. 

2.2.5 Thermodynamic Equations 

In this subsection, the thermodynamic equations which govern the large 
deflection (static or dynamic response) of the continuum in the presence of 
thermal and mechanical effects are presented since it is desired that the 
present general formulation include thermal effects for future use, although 
the illustrative example problems discussed herein do not include them. These 
thermodynamic equations will be used for the derivation of the material elastic 
and elastic-plastic constitutive relations presented in the next subsection. 

Consider that the continuum initially in the undeformed state with 
volume V q and bounded by surface A q is subjected to mechanical and/or thermal 
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loading. After deformation, its volume becomes V and is bounded by surface A. 


The first law of thermodynamics or the law of conservation of energy 
may be written as: 

For the Continuum as a Whole: 


[([ p„ ( F- ^ + B ) l /*ff vJA'-(( a ■ n„dA. =f(fp.0dV. (2 . 

A°<r Vo 


53 ) 


At any Generic Point of the Continuum : 

5*'% -Qu + P.B = p.U 

where B is the rate of heat input per unit mass 

Q is the rate of the heat flux vector across the body 
surface measured per unit undeformed area and 
Q * Q?g. m Q ^ 

A is the portion of the surface over which the heat 
flux proceeds 

U is the internal energy per unit mass 

p,F,T,A.n, S" j , and Y. . are defined as previously stated, 
o o oo o 13 

The second law of thermodynamics (or the increase of entropy 
principle) may be written as: 

For the Continuum as a Whole : 

/// ft s dv. -((( ft f <w. Q ^r dA ° ? 0 


( 2 . 54 ) 


Vo 


Vo 


ha\\ 


( 2 . 55 ) 


At any Generic Point of the Continuum 

p„S-P.T + t 2 


( 2 . 56 ) 
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where S is the entropy per unit mass 
T is the absolute temperature 

Let the Helmoltz function H be defined as 


H = U- T5 

• Then, substituting for U of Eq. 2.57 into Eq. 2.54, results in 

S xi V*i -Qii+P'B = P„ (H + i'S' t 5T) 


(2.57) 


(2.58) 


It is convenient to eliminate p B from Eqs. 2.56 and 2.58 to obtain the in- 

o 

equality . ( , 

5‘%- -mh+ts)- ? q T„ ?o 

Alternatively, the thermodynamic equations can be expressed in a form 
which refers to the deformed state of the continuum. 


These thermodynamic equations, if desired, can be employed to deduce 
the equation of equilibrium, the continuity equation, and the boundary condi- 
tions. 


It should be noted that all of the equations derived so far are exact for 
either finite strain or infinitesimal strain; that is, no restriction on the 
magnitude of strain has been made. 


2.3 Constitutive Relations 

This subsection is concerned with the elastic and elastic-plastic materi- 
al behavior of the continuum; also, thermal effects and strain-rate effects 
will be discussed. 

The elasticity theory (linear or nonlinear) postulates that a unique 
relation exists between the instantaneous states of stress and strain. When 
the material is deformed into the plastic range, a unique relation does not 
hold in general between stress and strain, but a functional relation exists: 
the strain depends not only on the current state of stress but also on the 
history of loading. 

At finite strain and/or elevated temperatures , the coupling influences 
occurring through thermal effects, elastic deformation, and plastic flow may be 
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significant in some cases, and an initially isotropic and homogeneous material 
may exhibit appreciably anisotropic and nonuniform behavior. The coupling in- 
fluences include for example: (1) the elastic modulii of the material will 
change as functions of temperature and continued plastic straining; (2) heat, 
which is generated from the dissipation of plastic work can induce the redis- 
tribution of temperature; also, (3) the Bauschinger effect and hysteresis loop 
are observed which are believed to be due to differential hardening of the 
variously-oriented crystals. These complex coupling functional relations, and 
the limited available experimental data make it very difficult for one to carry 
out the finite strain and/or high temperature elastic-plastic analysis with con- 
fidence. 

In classical thermal-elastic-plastic theory (Re£s. 13, 14) which is re- 
stricted to infinitesimal strain, considerable simplification is attained by 
decoupling the effect of thermal elasticity, mechanical elastic deformation, 
and mechanical plastic flow. The total strain is assumed to be decomposed 
into elastic (thermodynamically reversible) strain, and plastic (thermodynami- 
cally irreversible) strain. The elastic strain can further be decomposed into 
thermal and mechanical parts. Each component has the same invariant properties 
as the total strain, but is not expressible in terms of the displacements. Only 
the total strain can be related to the displacements. The elastic strain is 
related to the stress through a linear- function relation, and the thermal- 
elastic strain is treated as an initial (or a prescribed) strain. The elastic 
modulus and heat capacity are considered to be invariant under plastic flow; 
that is, the influence of plastic flow on the thermal-elastic characteristics 
is neglected. Only very recently, some pertinent work has been extended to 
finite strain (Refs. 15, 16). However, these theories do not completely 
agree with each other, and further research is required. 

In order to make the present study self-contained, the thermoelastic 
constitutive relation is presented in Subsection 2.3.1. Thermal-elastic- 
plastic material behavior is discussed in Subsection 2.3.2. Subsection 2.3.3 
is concerned with the strain-rate effect. The general derivations are based 
on the thermodynamic equations and follow those of Refs. 14 and 15. 
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2.3.1 Thermoelastic Material Behavior 

Consider an elastic continuum initially in a reference state with neither 

stress nor strain, and of uniform absolute temperature, T q . This continuum is 

then subjected to a temperature distribution T(^ 1 ,t) as well as body forces and 

surface tractions which may vary both spatially and temporally. Now let it be 

assumed that the Helmholtz function H, entropy S, heat flux Q and stress S 1 " 1 

K 

are functions of strain y. . , temperature T, and/or temperature gradient T 

i j ,k 

(Ref. 15). Thus, 

H = H C 'by, T ) 

S = 5 Chy - T ) 

=i W (Y, y,T) 

* w 

Q k = T, Tj ) 

(2.60) 


H = 5 =5=0 


when = 0 and "J*= “fo 


(2.60a) 


Then the thermodynamic equation, Eq. 2.59, becomes 




(2.61) 


At a given state of Yj,^ and T, the quantities y.^ and T can be chosen arbi- 
trarily. Hence, it may be concluded that 


5 = Po 


o iJd. 

0 dj 


Q T. k 


(2.62) 


If the Helmholtz function, H, is specified (for example in terms of 
and T) , the general thermal-elastic relation may be derived from Eq. 2.62. 

For infinitesimal strain, the elastic relation between stress and strain 
may be assumed to be linear and is known as Hooke's law. Thus, the Helmholtz 
function H may be expressed as a quadratic function of Yj.^ (cubic and higher 
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(2.63) 


order terms are neglected) : 


H=re,‘ >u Vi» +E /%- + e s + o(V 

<L i k£ x "i 

where E^ , J , E^ are functions of T and 

E E /^ ^ V 7/c E ^ = £ ^ 


Then, from Eg. 2.62 one may obtain 

s«-E* , r u + £*‘ t 


(2.64) 


Let the tnermal strain be denoted by y T . ; since it is required that = 0 
^ a 3 

when y ■ • = Y , . , one has 

P *1 = _ C ^ T 

t z ~ 'fcj (2.65) 


and Eg. 2.64 becomes 


:" f ' = E/ ,a (V k< -'f w ) 


(2.66) 


Equation 2.66 gives the linear thermoelastic stress-strain relation. The 
thermal strain is treated as an "initial strain" or a prescribed strain which 
produces no stress. When the material is isotropic both elastically and 
thermally, one may choose 


-y k , = V $u 


(2.67) 


and (Ref. 17) 


E '' a =/I [G‘V' + 


( 2 . 68 ) 


where h , X are temperature-dependent Lame 7 constants. Further, if the relation 
between the thermal strain and temperature is assumed to be 


it follows that 


dy = d(T) d T 

k 

u 

y = f T d(T) dj 

JT n 


(2.69) 


(2.69a) 
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where a is the temperature-dependent thermal expansion coefficient. Then, for 
an isotropic and homogeneous material, the linear thermoelastic stress-strain 
relation can be written in the form 


5: 


_E_ 




I-IV 


si) 


E_ ri 
\-V> i 


1 


d d T 


To 


(2.70) 


wnere V is Poisson's ratio, E is Young's modulus, and y. is defined to be 

. mi 


T, 




(2.70a) 


It should be noted that the coupling between thermoelasticity and heat 
conduction is neglected in the above equations. 


2.3.2 Thermoelastic-Plastic Material Behavior 

There are two types of common plasticity theories, termed "flow" and 
"deformation". The deformation theory of plasticity assumes that, as in elas- 
ticity, there exists a one-to-one correspondence between stress and strain. 

The flow theory of plasticity states that there is a functional relation be- 
tween the incremental stress and the incremental strain. Only for propor- 
tional loading where the stress ratio remains constant, and for a certain 
restricted range of loading paths other than proportional loading (Ref. 18) 
(through the assumption of the possibility of singularity in the yield surface) 
does the deformation theory agree with the flow theory. In order to include 
the capability to analyze general loading paths including loading, unloading, 
and cyclic loading, the "flow-type" theory will be incorporated into the present 
analysis . 

The behavior of a general elastic-plastic material can be characterized 
by the following two ingredients. First, assume the existence of a boundary 
(yielding surface) in stress space which defines the elastic domain. Within 
the boundary the continuum deforms elastically. Only at the boundary, the 
onset of plastic flow (irreversible deformation in a thermodynamic sense) is 
possible and no meaning is associated with the region that is beyond the bound- 
ary. Second, one employs a flow rule which describes the behavior of the ma- 
terial after yielding has started; it gives the relation of plastic flow 
(strain increments) to the stress (or stress increment) and the loading history. 
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Another basic assumption in the theory of an elastic-plastic continuum 

is the introduction of a plastic strain tensor. The plastic strain tensor, 

yjj\ , is assumed to have the same invariance properties as does the Lagrangian 

strain tensor, y. .. The quantity y~ . is related to r , by an elastic strain 
e 3-3 13 13 

tensor y. ., in the form (Ref. 15): 

Xf =Xy -"Xxj (2 . 71) 

No kinematic meaning is given to yj\ and y? j , but only their sum, y^_. as de- 
fined by Eq. 2.8 is related to displacement by Eq. 2.9. 

The yield surface, 9, is assumed to be expressible in terms of certain 
variables and may be expressed as 

7 


$ (s‘\ , T, <) =0 


(2.72) 


where K is a hardening parameter which depends on the strain history. Differ- 
entiating Eq. 2.72 with respect to time gives 

7 


i, .i! ; 
§ - d ji S 
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dT 


? K 




(2.73) 


If the conditions (a) 9 <_ 0 and 9 = 0 or (b) 9 < 0 are satisfied, the state 
change can only be elastic (reversible) , any plastic deformation (which may 
have been incurred earlier) remains unchanged. Thus, 

■7 d $ d 3? • _ 

^ when ^ + “J"<0 and = 0 (unloading) 


or 


and $=0 (neutral loading) 


§ < 0 


(elastic deformation) (2.74) 


If it is postulated that the plastic strain rate y?^ is linearly related to 

S 13 and T and with the consideration that y? . is zero for a neutral change of 

13 

loading, the following linear relation may be chosen as a reasonable approxima- 
tion: 
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X, = MtF' 


(2.75) 


dL cw + Li-f 
?5 wS n T?0 


$ = 0 


(loading) 


where D. . is a symmetric function of S k *’ f T, and the previous strain history, 

13 • kJl • 

but not of S and T. This assumption, which is suggested by Ref. 14, is 
based on the consideration that in a crystal grain, when the stress state is 
such that the resolved shears along certain slip direction reach the critical 
value, the plastic-strain increment occurs. Hence, as a statistical average 
over all grains, a definite macroscopic stress is needed. The stress-increment 
only determines the magnitude of the plastic-strain increment. Further, the 
(thermodynamic) Helmholtz function H, entropy S, and heat flux Q^, may be 
assumed to be functions of y® , yp , and T (Ref. 15). Thus, 


= H CY W 

= $ 

= q'o; 


, T) 

-ril, T, T,k) 


(2.76) 


Then, from the thermodynamic equation, Eq. 2.59, one may deduce that 


Q d H 


(2.77) 


If the functions $, D „ , and H are Chosen, Eqs. 2.74 and 2.75 will per- 
mit the determination of the plastic strain increment corresponding to the 
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stress increment and temperature increment, and Eg. 2.77 predicts the relation 
between elastic strain components and the stresses. Zt should be noted that 
no restriction is made on the magnitude of the strain and the symmetry of 
material properties in the above equations. But, for finite strain the assump- 
tion that the total strain can be linearly decomposed into elastic and plastic 
components and the postulate that plastic strain rate is linearly related to 
stress rate and/or temperature rate requires careful scrutiny, as claimed by 
Ref. 16, because of the non-Euclidean characteristic of the intermediate 
stress-free state associated with the plastic flow. However, for infinitesimal 
strain, these assumptions constitute a satisfactory approximation, since the 
deformation behavior exhibited by most engineering materials in infinitesimal 
strain is that the irreversible plastic flow causes no volume change and the 
reversible elastic strain is small. The specific considerations upon which the 
functions $ and D_ are chosen will be discussed in the following. 

In simple tension and/or compression tests of most engineering materials, 

the uniaxial stress-strain relation indicates the existence of a yield stress 

(elastic limit) , a , beyond which plastic (permanent) deformation takes place, 
o 

However, under multiaxial stress (and strain) the behavior is much more compli- 
cated. Various yield criteria and flow rules have been proposed for the pre- 
diction of the onset of plastic flow and the relation among plastic flow, stress, 
and stress history. Among them is the Mises-Hencky yield criterion and its as- 
sociated flow rule which usually fits experimental observations better than the 
Tresca criterion, for instance, for polycrystalline metals and yet is mathema- 
tically simple. The Mises-Hencky rules will be discussed and adopted in the 
present analysis. 

The Mises-Hencky yield criterion may be physically interpreted as 
"yielding begins whenever the distortion energy equals the distortion energy 
at yield in simple tension". Thus, hydrostatic pressure, for an isotropic 
material, (tension or compression) does not affect the yielding, plastic flow, 
and resultant hardening. Stated otherwise, no plastic work is done by the 
hydrostatic component of the applied stress. This implies that there is no 
plastic (or irreversible) change in volume. Thus, 
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-r* =° 


*4 - 0 


(2.78) 


* sp 

where is the spherical plastic strain-rate tensor. . For an initially 
isotropic material, the Mises-Hencky yield function can be written in the form 

(2.79) 

X> . u . I i ^ I .2 

2 i x 


where 


# = l t - jcrJ<T) = o 


? _ I 4-i - 1 r T x T t 

I. - p ? I .3 1 J 


D. 

T. 1 

3 


is the deviatoric stress tensor of T. in mixed 

3 

component form 


CT q (T) is the temperature-dependent yield stress in 
simple tension. 


(2.79a) 


Equation 2.79 represents a hypersurface in nine-dimensional stress space. Any 
point on this surface represents a point at which yield can begin. Also, it 
is assumed that the temperature effect changes only the size of the yield sur- 
face (through O q ( T)) and that the surface will retain the same shape and orienta- 
tion as in the initial reference temperature state. For moderate hydrostatic 
pressure and temperature, this is a very close approximation and is found to be 
in good agreement with experiment for common metals. 


Then, from Eq. 2.78, the function in Eq. 2.75 must satisfy the 
following two necessary conditions: (1) the restriction = 0 must be im- 
posed to ensure zero plastic volume change and (2) the principal axes of the 
plastic strain increment and of D„ must coincide with the principal axes of 
stress, since the material is isotropic. These two conditions can be satisfied 
with sufficient generality by choosing 


= a 


lE 
dS ‘ 


(2.80) 


where A and S' are functions of the deviatoric stress invariant and possibly 
also of the stress history. 


Next, the derivation of the function 0 will be based on Drucker’s postu- 
late (Ref. 19). Drucker 1 s postulate states that "the net work performed by 
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an external agency over a cycle of loading and unloading is nonnegative, posi- 
tive if the plastic flow has occurred, and zero if there is only change in 
elastic strain”. This postulate leads to the following two requirements on the 
plastic behavior: (1) the convexity of the yield surface with respect to the 
origin of the stress space, and (2} the outward normality of the plastic strain 
increment vector to the yield surface. It is noted that the Mises-Hencky 
yield function given by Eq. 2.79 satisfies the convexity condition imposed on 
the yield surface. The second requirement suggests that the function D may 
take the same form as $; thus. 



art 


r iM 

l ? jkt 


S U + 



where 


A '=*($* 


and 


= A 


' a_? 

J~) is a nonnegative scalar quantity 


’■ is the deviatoric plastic strain 

t 

rate tensor. 


(2.81) 


Substituting the Mises-Hencky yield function, Eq. 2.79, into Eq. 2.81 and 
using Eq. 2.28, one obtains the following flow rule: 

■ k . cr. d ■ 


yii. f. 


(2.82) 


^ ■ 


If . 

si- 


where 


X ™ A(g/G) 

D. t 

S . is the deviatoric stress tensor of S . in 
3 3 

mixed component form and is defined to be 

s 1 =■ S 1 r i s“ <5* 

3 3 3 m 3 


(2.82a) 
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For certain materials, the yield surface will change in case of continued 
straining beyond the initial yield. The change of the yield surface (called 
subsequent yield function) that characterizes the work hardening (or strain 
hardening) behavior of the material depends on the loading history. There are 
several hardening rules available to describe the subsequent yield function. 
Among them, the popular ones are "isotropic hardening" and "kinematic hardening". 


Isotropic hardening assumes that during subsequent yielding from a plas- 
tic state, the yield surface will expand uniformly with respect to the origin 
in the stress space but will retain the same shape and orientation as it had 
initially. It does not take into account the Bauschinger effect that is ex- 
hibited by many structural materials. Mathematically, the subsequent yield 
function for an isotropic hardening material can be put in the form: 


T D -2 

i- I, - ? <r = ° 


and 


cr = <r 


fY W, T) 


(2.83) 


where dw 1 ’ « dy^_. > 0 is the irreversible plastic work expended. 

To account for the Bauschinger effect, Prager introduced the "kinematic 
hardening rule" which postulates that during subsequent plastic flow, the yield 
surface translates (as a rigid body) in the stress space and that it will retain 
the same size, shape, and orientation that it had initially. Mathematically, 
this caui be expressed as 

si x 'i 


$ ( T - <0=0 


(2.84) 


where (y^_. ) which represents the translation of the referenced origin 

in the stress space of the yield surface and depends on the degree of hardening . 


It should be noted that, in the present analysis, the strain-hardening 
behavior of the material will be accounted for by using the well-known 
"mechanical sublayer model” (Refs. 20, 21, 22), The feature of this model is the 
inclusion of kinematic hardening (see Appendix A). In this model, the material 
at any point is conceived of as consisting of "sublayers"; each sublayer behaves 
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as an elastic, perfectly-plastic medium, having common strain and a common 
elastic modulus but appropriately different yield stresses. For a perfectly- 
plastic medium (nonhardening medium) , the yield surface remains unchanged in 
the case of continued straining beyond initial yield. 

2.3.3 Strain-Rate Effect 

For some engineering materials, the uniaxial stress-strain curve of 
dynamic loading tests will usually be different from those of static tests. 
The yield stress frequently increases with the increase of strain rate. How- 
ever, this strain- rate effect arises in a less understood fashion. 

The most general form of the uniaxial stress-strain relation including 
the strain-rate effect may be expressed as (Ref. 23) 

£ = X ( C, Z, T) + Y (O’, £, f ) (2.85! 

where 

a and e are the uniaxial stress and strain, respectively 
T is temperature 
X and Y are functions of state 


Under static loading tests, the rate-term may be neglected and one has the 
static uniaxial stress-strain relation 


Y ( <r, e , T ) = o 


( 2 . 86 ) 


One of the simple methods for approximating the strain-rate effect and 
which is in good agreement with experiment on certain types of common metals 
and alloys is to assume that the uniaxial stress-strain curve is affected by 
strain-rate only by a quasi-steady increase in the yield stress above the "static" 
test yield stress and that the elastic deformation is independent of strain rate. 
The increase in the yield stress under strain rate may be expressed in the follow- 
ing simple form (Ref. 24) : 



(2.87) 
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where 


a ' is the static yield stress 

£ is the uniaxial (or equivalent) strain rate 
D and p are empirically-determined constants for the 


a 

y 


material, which may be temperature dependent 
is the yield stress under £ 


Alternate approximations are discussed, for example , in Refs . 25 and 26. 
Equation 2.87, however, has been employed with satisfaction in recent years, 
and is adopted for the purposes of the present study. 
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SECTION 3 


GENERAL OUTLINE OF METHOD OF ANALYSIS 

3.1 Background and Survey of Pertinent Literature 

In this investigation, attention is restricted to methods for analyzing 
dynamic structural response, with principal attention devoted to the transient 
responses of structures which are subjected to transient external loads such 
as those arising from gusts, blast, impact, etc.; simple-harmonic vibratory 
response is of secondary interest. Explicitly excluded from consideration is 
the "short time" or "early time" response which is often called "material re- 
sponse", and which pertains to the nature, propagation, and effects of intense 
stress waves in the material as a result of severe impact or impulsive loads 
(mechanical and/or thermomechanical) applied to the structure; roughly the 
time span of interest for this type of response is of the order of from 1 to 
100 microseconds. Only the "late time” response which is usually termed 
"structural response" (in contrast with “material response”) is discussed here; 
such responses involve times of interest extending from time zero to 1 milli- 
second or perhaps to several hundred milliseconds; this type of response per- 
tains to the transient bending and/or stretching behavior of overall structures 
or of structural components such as beams, rings, plates, shell panels, etc. 
Furthermore, principal interest in this study centers upon transient structural 
responses involving both large deformations and elastic-plastic material be- 
havior. Sought is information on both the peak transient responses (deflections, 
strains, and stresses at any selected point in the structure) together with the 
time of occurrence of that peak and the permanent deformation condition of the 
structure after subsidence of the externally-applied transient loading. 

Accordingly, it is perhaps useful to review briefly the available analy- 
sis methods for various categories of transient structural response problems . 
Convenient categories are indicated in the following: 
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Further, it is useful to tabulate, as follows, the principal available transient 

< 4 > 4 ‘ 

response analysis methods and to indicate which are applicable to each of the 
problem categories cited above. 


Transient Response 
Analysis Method* 

Applicable Structural 
Response Category 

I II III IV 

(SD-LE) (LD-LE) (SD-EP) (LD-EP) 

A. Normal Mode Method (NMM) 


(Refs. 27-29, for example) 

X X 

B. Assumed Mode Method (AW) 


1. (Refs. 27, 28) 

X 

2. (Refs. 30,31) 

X 

3. (Ref. 31) 

X 

4. (Ref. 31) 

X 

C. Conventional humped Parameter (CLP) 


1. (Refs. 27, 28, 32) 

X 

/ 2 . (Refs. 33, 34) 

X 

3. (Refs. 33-35) 

X 

4. (Refs. 33-35) 

X 


+ This category includes strain hardening as well as strain-rate dependent and 
temperature-dependent material behavior. 

Some special-purpose less comprehensive methods are included in a separate 

brief discussion later. 

* 

Only typical references in each category are cited; the citation of a complete 
bibliography is not intended. Also, note that there are many uncited references 
in each category, for which only static responses are discussed; selected such 
references are discussed at pertinent places later in the text. 
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Transient Response 
Analysis Method 

Applicable Structural 
Response Category 

I II III IV 

(SD-LE) (LD-LE) (SD-EP) (LD-EP) 

D. Finite Difference in Space and Time (FD) 


1. (Refs. 2, 3, 36-41) 

X 

2. (Ref. 42) 

X 

3. (Same as under D.4) 

X 

4. (Refs. 3, 43-60) 

X 

E. Forward Integration in Space 


with Finite Differences in Time 


1. (Refs. 4, 5, 61-62) 

X 

2. (Ref. 63 + ) 

X 

3. (No Refs.) 

** 

4. (No Refs.) 

** 

F. Finite Element in Space 


with Finite Difference in Time 


(FES-FDT) 


1. (Refs. 64-73) 

X 

2. (Refs. 74-76) 

X 

3. (Same as under F.4) 

X 

4. (Refs. 77-78) 

X 

G. Finite Element in both 


Space and Time (FES-FET) 


1. (Refs. 12, 79) 

x ** ** ** 


In addition one could have various combinations of the above. For example, one 
often hears speculation (Kefs. 80 and 81, for example) concerning combining 
finite difference and finite element procedures in space to take maximum ad- 
vantage of the special merits of each method for appropriate parts of the struc- 
ture (that is, use finite differences for smoothly-varying regions of the struc- 
ture, and finite elements in regions of structural irregularities such as ir- 
regular cutouts, etc.) — in combination with an appropriate finite-difference 
_ 

Method applicable but not described explicitly in the literature. 

+ Only static problems are discussed. 
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operator in time. However, a concrete illustration of this approach in the 
open literature has not been found. 

For small-displacement, linear-elastic structural behavior (Category I) , 
the governing equations are linear and hence the use of the normal mode ap- 
proach together with solution superposition is an applicable and efficient 
method of solution. Each of the other cited methods (B through G) may also be 
employed for this problem category but these ere often more laborious and less 
attractive . 

When large displacements are encountered but the material behavior is 
still linear-elastic (Category IZ) , the governing equations become nonlinear. 

In this case the "normal modes" become amplitude-dependent; thus, the normal 
mode method of analysis becomes impractical even if one approximates the be- 
havior as being piecewise linear. For such cases, methods B through G may be 
used. 

On the other hand, for problems involving small-displacement, elastic- 
plastic behavior (Category III) , one can legitimately regard the system as re- 
taining its original small-displacement linear-elastic identity insofar as 
its normal modes are concerned (as identified by the mass and stiffness charac- 
teristics of the structure) . Plasticity effects may be taken into account as 
equivalent plastic forces which are amplitude dependent. Hence, the NMM becomes 
unwieldy since "solution superposition" is destroyed by the essential nonlinear 
character of this type of problem. Accordingly, methods B, C, D, F, and G are 
better suited for this Category III type of problem, with methods D and F most 
widely employed. 

Finally, Category IV poses the most severely nonlinear set of conditions 
both geometric nonlinearities and material nonlinearities are present. While 
methods B, C, D, F, and G are all applicable to this type of problem, only the 
finite-difference method (FD: method D) has been highly developed for analyzing 
this type of transient response problem (Refs. 3 and 43-60) — only Refs. 77 
and 78 have reported the use of the finite-element method (FE; method F) for 
this transient-response-problem category. 

The finite-difference method is a very general and powerful method which 
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has undergone more extensive development for each category (I, II, III, and 
IV) of problems than any of the other methods mentioned — much literature, 
not cited here, exists for both static and transient response problems. In the 
past 15 years , however , the finite-element method has undergone much develop- 
ment for both static and transient response problems , especially those in 
Category I, with successively lesser developments in categories II, III, and IV. 

For the finite-element analysis of static problems in Category II (large- 
deflection, linear-elastic) , Turner et al. (Ref. 82) was the first to report, 
and many extensions have since been reported. Among these, four general schemes 
have emerged: 

1. The first approach is based on establishing a linearized incre- 
mental formulation; the load is applied incrementally. This 
geometrically nonlinear problem is solved by carrying out a 
re-evaluation of the element-stiffness matrix at each load 
stage. This "incremental change" in the stiffness matrix is 
termed the "geometric stiffness matrix" (Refs. 83-88) . 

2. Schmit et al. (Ref. 89) has solved the geometrically nonlinear 
problem by seeking the minimum of the total potential energy 
function by a direct energy search procedure. 

3. In the third approach, the governing equations are solved by 
an iteration procedure such as the Newton-Raphson technique 
(Refs. 90, 91). 

4. Stricklin et al. (Refs. 92, 93) analyze this geometrically non- 
linear problem by treating the large deflection terms as equiva- 
lent force terms which are derived from the pertinent energy ex- 
pression in the variational formulation used; for those special 
terms, a restricted assumed displacement field is used for the 
finite-element analysis to avoid certain numerical difficulties. 

This approach has also been extended to transient problems (see 
Refs. 75-76) ; when an implicit finite-difference time operator 
is used, iterative and/or approximating extrapolative procedures 
are needed to take those terms into account. 
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For the finite-element analysis of static problems in category III 
(small-deflection, elastic-plastic) , various investigations have been carried 
out (for example. Refs. 94-99). Two means for treating material elastic-plastic 
behavior were employed: 

1. The method of initial strain (see Mendelson, Ref. 100). 

The idea is to modify the load-deflection equations of 
equilibrium so that the linear-elastic stiffness matrix 
is employed throughout the elastic-plastic load-deflec- 
tion range of interest, and plastic effects are taken 
into account through the use of "effective plastic loading" 

(Refs. 94-97). 

2. The tangent modulus method. This method is based upon the 
linearity of the incremental laws of plasticity, and deals 
with this behavior in a piecewise linear fashion (Refs. 98,99). 

The relative merits of these two approaches have been discussed by Marcal 
(Ref. 101) . 

In Refs. 102 and 103, for example, the finite-element method is discussed 
for static problems in Category IV (both large-deflections and elastic-plastic 
behavior). Also in Ref. 104, the finite-element equations of dynamic equilibrium 
and an increment stiffness equation are derived; however, no transient response 
predictions are reported. Reference 102 utilizes the initial strain concept 
for including plastic behavior, while Refs. 103 and 104 employ the tangent 
modulus approach; in all of these cases, a linearized increment formulation re- 
sults. 

Recently, Salus, Ip, and VanDer linden (Ref. 77) have described the formu- 
lation and application of a finite-element approach for predicting the large- 
deflection elastic-plastic transient response behavior of beam-type structures. 

The resulting predictions, for several examples reported, are in good agreement 
with pertinent finite-difference predictions and experimental results. This 
formulation is of the assumed displacement type but is not based upon variational 
principles. Emphasizing the plastic part of the behavior, only linear displacement 
fields are introduced for all of the displacements , and transverse shear 


38 



deformation is included. While elastic, perfectly-plastic behavior is taken 
into account, the effect of transverse shear on yielding and flow appears to 
be accounted for only in an indirect fashion. 

From the preceding brief review of the finite-element literature, it 
is apparent that only limited finite-element developments have been reported 
for Category XV transient response problems: involving both geometric non- 
linearities (large deflections) and material nonlinearities (elastic-plastic 
behavior). Accordingly, the present study is devoted to developing, evaluating, 
and applying a finite-element variational method for analyzing the large-deflec- 
tion elastic-plastic transient and permanent deformations of simple structures. 

Finally, note should be made of a group of simpler but more restrictive 
methods for estimating transient responses and/or permanent deformations of 
structures; these include the "rigid-plastic" method and the "energy-absorption" 
method. In the rigid-plastic method (Refs. 105-116), the elastic deformation 
behavior is assumed to be negligibly small compared with that arising because 
of plasticity. Using the concept of plastic hinges, the equations of motion 
for a structure are derived for various levels of approximation for the kine- 
matics of the system; the deflections may be small or large and the material 
may be rigid-plastic, rigid-visco-plastic, and/or strain-rate dependent rigid- 
plastic. Except for small displacements, the governing equations are nonlinear 
and numerical solution methods must be used. Often with this method one ob- 
tains reasonably good estimates of the permanent deformation; however, the 
transient response is always rather badly in error: the peak deformation is 
always under-estimated and the time to peak response is invariably too great. 

For "trends and/or parametric studies" wherein trends rather than accurate 
values are desired, the rigid-plastic approach is often useful. Various 
"bound theorems" have been developed (Refs. 117-122) to estimate upper bounds 
and lower bounds on the deformations of impulsively loaded rigid-plastic and 
elastic-plastic structures; these form useful supplements to the rigid-plastic 
analyses . 

The "energy absorption" method (Refs. 123-127) is used only to estimate 
the permanent deformation of the structure; no transient response information 
is obtained. In the energy-absorption approach, it is assumed that the primary 
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type or pattern of inelastic deformation is known for the structure being 
analyzed. Also, it is often assumed that the elastic energy stored in the de- 
formed structure is negligible compared with that absorbed by plastic work. 
Accordingly, the permanent deformation is estimated by equating the plastic 
work absorbed by the structure in deforming in the prescribed deformation 
pattern to the total work done on the structure by the externally-applied loads. 
Thus, for example, if a structure were loaded impulsively, the structure would 
be given an initial kinetic energy which will, in general, would produce both 
rigid-body motion and structural deformations. The portion of this kinetic 
energy which is available to produce structural deformations is then equated to 
the plastic work absorbed by the structure — thus providing an estimate of the 
permanent deformation state of the structure. The success of this method de- 
pends upon reasonable a priori estimates of the primary pattern of deformation 
occurring during the large-deformation response; for complex structures sub- 
jected to arbitrary transient external loading, this a priori knowledge is 
usually lacking. For certain structural configurations and external loadings, 
however, accumulated experimental evidence is available for making such de- 
formation-pattern estimates. In such cases, very reasonable permanent deforma- 
tion estimates result; Greenspon (Refs. 124-127) has made very effective use 
of this excellent and efficient special-purpose approach. 

The deficiencies of the rigid-plastic method and the energy absorption 
method have been overcome by employing numerical methods such as the finite- 
difference (FDS-FDT) method and the finite-element (FES-FDT) method. The 
present investigation has as one of the main objectives, the extending of the 
FES-FDT approach to analyze in an accurate and rigorous manner the large- 
deformation elastic-plastic transient and permanent deformation of transiently 
loaded simple structures. The intent is that, if desirable, the techniques 
developed in this study for simple structures could be extended to analyze 
more complex structures. 

3.2 Variational Derivation of the Equations of Motion 
Utilizing the Finite-Element Approximation 

As noted earlier, the finite-element approach is utilized in conjunction 
with the Principle of Virtual Work and D'Alembert's Principle (in short: 
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PVW-DAP) integrated over time + between times t^ and t as in Hamilton's Principle 
to obtain the equations of motion for a continuum or structure which is permitted 
to undergo large-deflection elastic-plastic transient deformations. In the 
present investigation, the assumed displacement form of the finite-element method 
is employed; when this is applied via PVW-DAP and the usual reductions are carried 
out, the resulting equations of motion for an undamped system appear in the follow- 
ing conventional form (Refs. 93, 94, 96, 128); 


=!FVfF;Vr !) i 




if;) 


(3.1) 


where 


q* (q*) 

[M] 

[K] 

{F*} 

{F^(q* 2 ,q* 3 )} 


(f* l > 


cure the global generalized displacements (accelera- 
tions) 

is the mass matrix for the complete assembled dis- 
cretized structure 

is the usual stiffness matrix (for linear elastic 
small displacement behavior) of the complete 
assembled discretized structure 
is the vector of externally-applied (global de- 
noted by a superscript *) loading 
represents a "generalized loads” vector arising 

from large deflections and is a function of 
2 3 

quadratic (q* ) and cubic (q* ) displacement 
terms — a nonlinear force contribution 
is a generalized loads vector arising from the 
presence of plastic and/or thermal strains, and 
is associated with the linear terms of the strain- 
displacement relations. 


The use of the time integrated variational statement is optional; the use of the 
PVW-DAP directly is sufficient since it holds at every instant of time. 
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(F^q*)} is a generalized loads vector of origin similar 

P 

to IF* } but is associated with the nonlinear 

p ' 

terms of the strain-displacement relations. 

Alternatively, by employing the finite element assumed-displacement ap- 
proach in conjunction with the FVW-DAP but by carrying out the reduction pro- 
cess differently, a computationally superior set of equations of motion is ob- 
tained. These "improved formulation" equations appear in the following form: 

where the quantities [M] , {q*}, {q*}, and {F*} retain the meanings given follow- 
ing Eq. 3.1. However, {p} can be shown to represent not only [K]{q*} of Eq. 3.1 
but also some plastic behavior contributions. Also the term [H]{q*} represents 
"generalized loads" arising from both large deflections and plastic (and/or 
thermal) strains. 

It is shown in Subsection 3.3.4 and in Section 5 that the improved formu- 
lation represented by Eq. 3.2 is much more efficient for a given solution accur- 
acy than is the conventional (Eq. 3.1) formulation. 

The improved formulation and the conventional formulation are developed 
in detail in Subsections 3.2.1 and 3.2.2, respectively in order to illustrate 
their similarities and differences. 

3.2.1 Improved Formulation 

In the finite-element-analysis method, the entire domain of the continuum 
is subdivided into a finite number of regions called "finite elements" or "dis- 
crete elements", each having a finite number of "nodes" as control points (see 
Fig. 3) . These nodes are usually located at the boundary of each element but 
may also be in the interior region of the element. The behavior of the actual 
continuum which has an infinite number of degrees of freedom is thereby de- 
scribed approximately in terms of a finite number of degrees of freedom at each 
of the finite number of nodes since the generalized displacements within each 
finite element is expressed in terms of (a) such variables called "generalized 
degrees of freedom" which are defined at the node points in conjunction with 
(b) suitably-selected interpolation functions to describe the distribution of 
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each quantity throughout the interior of each finite element. Applying this 
approach within the PVW-DAP framework results in a finite-sized system of 
second-order ordinary differential equations. The unknowns in these equations 
are the generalized degrees of freedom at each node of the complete assembled 
discretized structure (or continuum) . 

Although many different kinds of finite-element models exist (that is, 
displacement, equilibrium, hybrid, mixed, etc. — see Ref. 7, for example) , 
the assumed-displacement type of finite-element formulation or model has been 
chosen for development to analyze the present class of nonlinear transient 
response problems. A parallel study, of course, could be carried out by using 
perhaps each of the other types of finite-element models (see Appendix C) — 
and the relative merits of each could be assessed; such an undertaking, how- 
ever, is beyond the intended scope of this study. 

In the assumed-displacement- type of finite-element analysis , the gen- 
eralized displacements constitute the primary variables. Hence, one selects 
appropriate interpolation functions "anchored to” control-point values which 
are the nodal generalized displacements . In choosing appropriate interpola- 
tion functions for each finite element to be used in the assembled finite- 
element array, one may take into account the following "sufficient conditions" 
which will insure that the finite-element solution will converge to the exact 
solution as the continuum is more and more finely subdivided [129, 130] : 

1. Rigid-body modes of an element must be included; otherwise, 
the equilibrium conditions of the element as a whole will 
be violated. 

2. Uniform strain states must be included; otherwise, it cannot 
be assured that, as the mesh size is made finer, the strain 
will converge to the true state of deformation. 

3. The admissible conditions of compatibility should be satis- 
fied along the interelement boundaries as well as within 
the element. 

Accordingly, these guidelines are followed in selecting interpolation 
functions to represent continuous generalized displacements in the interior of 
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each finite element. A detailed description of this selection is given in 
Section 4 for an arbitrarily-curved beam element. 

Let it be assumed that the continuum or structure being analyzed has 
been subdivided conceptually into N finite elements. Thus, one may write the 
Principle of Virtual Work combined with P* Alembert's Principle integrated over 
time from t to t 2 , Eq. 2.52a, as the sum of the contributions from finite 
elements 1, 2, ... , N as follows : 



Jt = 0 


(3.3) 


where variations <5 are permitted only for the displacements, consistent with 
the displacement constraints for all times* except at and t^ and where for 
any element n: 


V„ 

K An 

(3.3b) 

SKn-fff dV " 

Vn 

(3.3c) 


+ At times t^ and t^, the configuration of the system is regarded as being known 
and hence the displacements at times t^ and t^ are not subject to variation. 
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In Eqs. 3.3 through 3.3c, V is the volume of the nth discrete element and A 

n n 

is the portion of the surface or boundary area of element n on which the surface 

traction is prescribed. Both V and A pertain to the initial undeformed 

n n 

configuration . The summation, E , extends over the "H" elements of the continuum. 
The other quantities have been defined following Eg. 2.52a. 

Let {v} represent the displacement field which consists of the components 
of the displacement vector expressed with respect to the undeformed base vector 
system of the coordinates (see Fig. 1). One chooses for each element an 
assumed displacement field of the form: 

M- 

where [U^ 3 )] is the matrix of appropriately assumed interpolation functions 
expressed in the coordinates ^ of a generic point within the element, and {6} 
represents a set of undetermined independent parameters which are functions of 
time only. 


It follows that the nodal generalized degrees of freedom which are the 
nodal generalized displacements, {q} , are defined in terms of the local coordi- 
nate system of each element and can be obtained by substituting the coordinates 
of the nodal points into Eq. 3.4. Accordingly, one may obtain 


f?l= [^J \P\ 


If one takes the same number of displacement parameters as the nodal generalized 
displacements, the transformation matrix [A] is a square matrix. By inverting 
Eq. 3.5, one has 


fP}= [Af'fji 


and Eq. 3.4 becomes 



where 


[ N] = [U<^](A]" 


Because IU(£ J )J and [A] are a priori chosen functions expressed in the 
coordinate only, they are not subject to variation; hence. 


fM- t N1 {Si\ 


Also, the time derivative of Eq. 3.7 becomes 


|vj = INj j 


(3.8) 


(3.9) 


By using Eqs. 2.9 and 3.7, one may obtain the corresponding strain 
at any point in the element as a function of position and the nodal generalized 
displacements, as follows: 


l t u A 1 

It follows that 


Syxf = l- D,j J li ^} + L V }■ *-■ D i J f' I 


(3.10) 


(3.11) 


where , and D_. are the appropriate associated differential operators 

which include both small and large deflection effects and which may be ex- 
pressed symbolically in the form* 

~ J L + 

l_D a i J = L 1 

l.D*J = L fsl% J 

t i (3.12) 

Employing Eqs. 3.8 through 3.12, Eq. 3.3 becomes 


* 

Explicit examples are given in Section 4. 
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n„i l + in nl^j s^v n \f\ 

Vn Vn 

-((I rNJ T p„/f)rfV„ -ff (N e ] T jJ^A.) 

- X [NJ T P, [NjrfV^jj^o 

K 


n = i 


(3.13) 


/ 


where subscript B is used to signify that the [N] are evaluated along the ele- 
ment boundaries. 

In more compact form. Eg. 3.13 may be written as 


'* , N 


; upi +[hitf \ - f f i ) 

J dt = o 


f\l 

Z 

n=- 1 


(3.14) 

where the following quantities are evaluated for each finite element* : 


lt\ = {({ fD J s‘*iV. 
v n 

[h]= fff ID^lD^j S‘ f dV n 
V n 

+ The evaluation of S 1 ^ is discussed in Subsection 3.3.2. 


(3.14a) 

(3.14b) 
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(3.14c) 


if j =/// [Nfp.ffi dv„ ■+([ rN 6 f|iH4 

V* a. 


[ m ) ~ f[f T NJ po [ Nj d 14, o. i4d) 

V, 

Since the element nodal generalized displacements {q} for different ele- 
ments are not completely independent, hence, a transformation is required to 
relate the element nodal displacements to a column of independent global 
(master) displacements for the discrete-element assembly, as follows: 


\h , b, b, - '".In) = 

(3.15) 


where are the element generalized coordinates in the local coordinate 
system, {q*} are the complete set of "master" generalized displacements in- 
cluding all node points of the complete assembled discretized structure, re- 
ferred to the global coordinate system. The quantity [J] includes the effect 
of transferring from local coordinates for each individual element to global 
reference coordinates for the system as a whole. When the coordinate trans- 
formation is not required (that is, when the element generalized displacements 
are already in the same direction as the global generalized displacements) , [J] 
is a simple Boolean matrix. 

Performing the indicated transformation, Eq. 3.14 can be written as 



Mfi ( f Pi + ) 

- LSfj — o 


(3.16) 
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where 


I Pi ' ( if it, ft. 


[H]= [J] 


hi 




IF*) = [ J] T fl', is, — • /» 1 
(MJ = (J] T 


W 2 



(3.16a) 

[J] 

(3.16b) 

i» 1 

(3.16c) 

] fJ] 

(3.16d) 


KM , 


Integrating the last term which appears under the integral sign of Eg. 3.16 by 
parts and using the fact that the virtual generalized displacements vanish at 
the two time limits, one obtains: 

)}dt = 0 ^ ^ 

Since the virtual global (master) generalized displacements {<5q*} are inde- 
pendent and arbitrary at each instant of time, the following equation of 
dynamic equilibrium results from Eq. 3.17: 


[mj \ t \ + (HjfrwF-j 


Given a set of initial conditions {q*}, {q*}, and {f*} at t = 0, and the proper 
boundary conditions, the system of second-order differential equations repre- 
sented by Eq. 3.18, may be solved in a step-by-step timewise fashion by using, 
for example, the finite-difference scheme. Further aspects of the solution 
process are noted in Subsection 3.3. 


Equation 3.18 represents the "improved formulation" form of the equa- 
tions of dynamic equilibrium. 


3.2.2 Conventional Formulation 

In order to make clear the source of the "apparent differences" between 
the improved formulation (Eqs. 3.2 or 3.18) and the conventional formulation 
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represented by Eq. 3.1, a parallel derivation of the latter is presented in 
the following. In fact, the derivations are identical up to and including 
Eqs. 3.13 and 3.14 through 3.14e. For the conventional assumed-displacement 
formulation, one proceeds to express the stresses in terms of the displace- 
ments via the stress-strain relations and the strain-displacement relations. 
That is, only the terms represented by Eqs. 3.14a and 3.14b undergo a further 
reduction by using the following stress-strain relation: 

C r "v "V ^ ) 

b “ t ( “ *fcl ' 13.19) 

where represents the components of the total plastic strain. Next, ex- 
pressing the total strain Y ki in terms of the generalized displacements {q} 
via Eq. 3.10, Eqs. 3.14a and 3.14b may be rewritten as: 


If) = Iff 

V„ 

i j ki 


if 

V„ 


= iki if} 


and (multiplying Eq. 3.14b by {q}) : 


(3.20a) 


V. 

- ffflD^LDfj E‘ /U < lD W J 
Vn 




(3.20b) 
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where 


CkJ =/([ ) j E i}k, LD kJJ dV„ 

V„ 

tf*} -gty) E‘ fi '( i <■ V M LD >^f\) W* 

Vn 

f^) = f 

p Vf, 


(3.20c) 


(3. 20d) 


(3.20e) 


f y' (3. 20f ) 


I f 7) =/f( [4 ,} iPJj e' , '“ ' rj dVJfi 

t 14, 


(3. 20g) 


and where superscripts L and NL refer to linear and nonlinear terms, respectively. 
Also note that the subscript ”p" has been applied to and to denote the 
presence especially of plastic strain but this can also refer to the presence 
of other "initial strains" such as thermal strains. 


L 


Substituting Egs. 3.20a and 3.20b into Eg. 3.14, one obtains 

( £ l ■ $ f j ( ( k ) l%\ - 1 f £ } - f ff } - / f p ■ - f f ) > 
£ lS£j 0]f ji J dt = o 


+ ' n=i 


n=i 


(3.21) 


where {f} and [m] are defined as previously stated, and 
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ff* L j 

{ ^ uv - 

* , 

-iPM e‘ ? M,j Jfj » I vv: l fl 

V„ 

(3.22) 


Then, transforming the element nodal displacements {q} to independent global 
displacements {q*} of the discrete-element assembly as described previously, 
Eq. 3.21 can be rewritten as 

f%xiKi it HFp-ifvHFthiH ■ )<* =0 <3. 2 3, 


where 


[K3 =[J1 T 


r kj , 

ki 


' is? / 


EJ3 


•#> 




/ F : l j=[j] T ii t L ,, ft,, ----.fM 


(3.23a) 


(3.23b) 


(3.23c) 


. _.+^i _ r _T » r ^ ^ 

'iff, if*, ■ i-pH j 


{ ri = t jfif' , i. 


CM] = CJJ 


rm. 


m 2 


Cn ; 


, f* } 
in 


(3.23d) 


(3.23e) 


(3. 23f } 
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Since the variations (_ 6q* J can be independent and arbitrary , the follow- 
ing conventional form of the equilibrium equations is obtained: 

[m jfrHKjffHF'WFp 

3.3 Timewise Solution of the Governing Equations 

In order to obtain the timewise solution of a set of equations of 
dynamic equilibrium such as Eq. 3.1 or Eq. 3.2, one may resort to analytical 
techniques or numerical techniques depending upon the mathematical (and/or 
physical) nature of the problem. 

For small-deflection linear-elastic behavior, for example, one may re- 
cast these equations into normal mode form and solve the resulting equations 
of motion analytically , mode by mode if the forcing functions are modally un- 
coupled or are properly sequentially coupled. Superposition of the forced 
response of each mode then provides the total response of the system. Alterna- 
tively, if desired, one may solve these equations (Eq. 3.1 or 3.2 in their 
present generalized coordinate form) timewise by using a finite-difference 
numerical procedure whereby one obtains a recurrence equation which provides 
a solution step-by-step in finite- time increments. 

Of these two methods, if the dynamic system is linear and is subjected 
to a transient forcing function which excites mainly the lower frequency modes 
of the system, it is frequently more convenient and efficient to use the 
normal mode approach. However, if the stiffness matrix (and/or the mass matrix) 
varies with time as in the present class of nonlinear problems, the normal modes 
also vary in time; hence, the normal mode approach becomes impractical. In such 
cases, the direct finite difference approach appears to be the only feasible 
method developed to date. Accordingly, the numerical finite-difference method 
is employed in the present study for solving equations of motion like Eq. 3.1 
or Eq. 3.2. 

In particular, the central-difference finite-difference time operator 
is employed for purposes of illustrating the solution process for the improved 
and for the conventional formulation in Subsections 3.3.1 and 3.3.2, respectively. 
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Then the comparative storage and computing requirements of these two formulations 
are discussed in Subsection 3.3.3. Next, since a wide variety of timewise finite 
difference operators or methods has been developed, a brief discussion of the ad- 
vantages and disadvantages of some of these methods is given in Subsection 3.3.4 
leading to the selection of the central-difference method for principal use in 
the present study. 

Finally, it should be noted that a somewhat different approach termed 
"the finite element method in space and time" has been developed (Befs. 12, 79). 
In this method, the initial value problem is handled as a boundary value prob- 
lem with the given initial conditions treated as prescribed boundary conditions. 
By using a procedure similar to that used in the discretization of only the 
space domain (as described in previous subsections) , the whole domain (both 
time and space) of interest is discretized into a finite number of discrete 
domains each having a finite number of nodes. At each node a finite number of 
degrees of freedom is permitted. Then, suitable interpolation functions, which 
are both time and space dependent, are selected throughout the interior of each 
discrete domain. By applying this approach within the framework of the time 
integrated PVW-D' Alembert variational statement results in a set of simultane- 
ous algebraic equations in the space and time degrees of freedom. Depending 
upon the type of interpolation function used (especially for the time domain) , 
one may obtain sets of equations corresponding to the use of many of the con- 
ventional finite-difference time operators, as well as a variety of other sets 
of simultaneous algebraic equations. In general, improved convergence and 
calculation stability will be observed as the entire space and/or time domain 
is more and more finely subdivided. 

3.3.1 Solution Process 

As indicated earlier, the equations of motion (Eqs. 3.1 or 3.18 and 
Eq. 3.2 or Eq. 3.24) are to be solved at a sequence of instants in time At 
apart by employing the following central-difference finite-difference approxi- 
mation for the acceleration "6 at any instant t : 
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% = 1 hiZls d -+ o (4t/ 

ym (4tr 


Also, one may approximate the velocity q_ at time t by 

to m 


(3.25a) 





2 (At) 


+ 0 ( A t ) 


(3.25b) 


Now note that at any time instant t , Eq. 3.24 for the conventional 


formulation may be written exactly as 


[M I f ft. - i K].f a - f Ft - f f;v Kl + 1 F T) 


-*ULi 


m (3.26) 


All quantities in Eq. 3.26 except [M) and (K) change with time (i.e., [MJ = 

= ^ = ^ = ••• ■ [K]^) . Assuming that all quantities in 

Eq. 3.26 except for {q*} are known at time t , one can solve Eq. 3.26 to 

m m 

obtain {q*} . Since one has already obtained the solution for {q*} at all 
m 

earlier time instants (that is, {q*} , {q*} ,, etc.), one can determine 

m m— i. 

approximately from Eq. 3.25a as: 




(3.27) 


where all quantities on the right-hand side of Eq. 3.27 are known. Thus, 

Eq. 3.27 gives an explicit evaluation for the generalized displacements at 

time t in terms of known information at times t and t Accordingly, 

m+1 m m-1 3 1 

Eqs. 3.25a and 3.25b are known as explicit finite-difference operators. 


Similarly, Eq. 3.24 for the improved formulation may be written at 

time t as 
m 
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(3.28) 


In Eg. 3.28, all quantities except [M] change, in general, with time. If the 

solution has been obtained for earlier instants in time , one may compute 

{q*} from Eq. 3.28 and then use Eq. 3.27 to obtain {q*} , . 

m m+1 

Assuming that at t = 0, the structure is at a known condition such as 

{q*} o = and {q*J o = {a}, one can readily obtain {q*}^ from 


(3.29) 


since {f*}^ is prescribed and all other quantities are known. 


In the timewise step-by-step solution procedure involving large-deflec- 
tion elastic-plastic transient responses via the spatial finite-element pro- 
cedure described in Subsection 3.2 (from which Eqs. 3.18 and 3.24 resulted), 
certain quantities in the governing equations change with time and hence must 
be re-evaluated, in general, at each instant in time. For example, for the 

improved formulation, Eq. 3.18 or 3.28, it is necessary at time t to evaluate 

m 

{p> m and [h]^ by using Eqs. 3.14a and 3.14b for each finite element — the 

"assembly" of this information then provides {p} and [H] via Eqs. 3.16a and 

in zn 

3.16b, respectively. Similarly, for the conventional formulation, one must 

evaluate {f HL } , {f L } , and from Eqs. 3.22, 3.20e, and 3.20g, respec- 

q m pm pm 

tively, for each finite element — the assembly of this information then pro- 


vides {f* NL } , {F*^} , and {F* wi '} via Eqs. 3.23b, 3.23c, and 3.23d, respec- 
qmpmpm 

tively. 

It is seen that the evaluation of {p} . [hj , {f* 1 } , {f L } , and {f NL } 

m' m q m pm pm 

involves volume integrals of certain quantities. For a structure undergoing 
large-deflection elastic-plastic behavior, it is impractical to evaluate these 
volume integrals analytically; instead, it is convenient and practical to per- 
form this integration numerically. Among the various numerical integration 
(or quadrature) methods, Gaussian quadrature (Ref. 131) appears to be the most 
efficient for a given accuracy. Accordingly, Gaussian quadrature is employed 
herein — this requires that the stresses and strains be evaluated at a 
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selected finite number of Gaussian stations or points over the "spanwise" and 
depthwise region of each finite element. 


Another ingredient which is common to both the improved formulation and 

the conventional formulation is the solution for {g*} at each instant of time 

zn 

t via Eq. 3.28 and Eq. 3.26, respectively. These equations are of the formi 


[mj {*">1 = l bit i 


for 


= 1 . £, 3 ,- 


(3.30) 


where 


(Ml 


{x(t)} 


m 


{b (t) } is 
m 


is a known banded positive definite symmetric 
matrix (the mass matrix for the restrained or 
unrestrained structure, whichever case is be- 
ing treated) 

is a vector of unknowns which must be de- 
termined by solving Eg. 3.30 
is a known vector (representing all terms 

except [M] {g*} in Eg. 3.26 or 3.28) 
sa 

-1 


In principle, one can always form the inverse matrix [M] and pre-multiply 
Eg. 3.30 by [M] to obtain 


[ M ) " [ M J = [Mf 

solution: 


which results in the solution: 


(3.31) 


since [M] ^ [M] = (II where [I] is the unit diagonal matrix. However, it has 
been found that independent of the number of time instants at which one wishes 
to solve Eg. 3.30 such a procedure is not as efficient as is the Choleski 
method (Ref. 132) . 


Briefly, the Choleski method involves factoring the matrix [M] to form 

a lower triangular matrix [L] and an upper triangular matrix (which is the 

T T 

transpose of the former) such that (M) = [L] [L] where [L] is the transpose 
of [L] . Thus, Eg. 3.30 may be rewritten as 

T 



(3.33) 


Next, form an intermediate matrix {y} which is defined as 

fy)„ - 

From Egs. 3.32 and 3.33, it follows that 

[L] I y}„= |b (t l 

At each time instant, one solves Eq. 3.34 for {y} very readily because [I>] is 

m 

a lower triangular matrix. One then solves Eq. 3.33 for {x} very rapidly also 

IQ 

by algebraic back-substitution. 

In Subsections 3. 3.1.1 and 3. 3. 1.2, the principal steps involved in the 
timewise solution of, respectively, (a) Eq. 3.28 for the improved formulation 
and (b) Eq. 3.26 for the conventional formulation through the use of the ex- 
plicit central -difference finite-difference operator (Eqs. 3.25a, 3.25b, and/or 
3.25c) are described briefly. This description is intended not only to point 
out the salient features of the solution but also to make clear the nature and 
extent of the similarities* and differences* between the improved and the 
conventional formulation. 


3 ■ 3 . 1 . 1 Improved Formulation 


The structure is represented by an assemblage of a finite number of 
discrete elements (also called finite elements) , and the geometric properties 
of each element are defined so as to approximate the actual geometry of the 
structure as closely as desired and/or feasible. The mechanical properties of 
the structural material are assumed to be known as a function of temperatures 
and strain rate. The structure is assumed to be subjected to externally- 
applied loads which are prescribed in both space (over the surface and/or 
throughout the volume of the structure) and time. The equations of motion, 

Eq. 3.18, to be solved according to the improved formulation are restated here 
for convenience: 


[M] if} - fpj -rHxri = if'} 


These aspects are then summarized in Subsection 3.3.3. 



Starting from a set of given initial conditions at time t = t = 0 on 

the generalized displacements ({q*} = {o}» for example} and the generalized 
• • * 

velocities {q} , one can solve Eq. 3.28 for {q*} at time t and then employ 
o o o 

Eq. 3.29 to compute {q*} ( . A slightly different but similar procedure is then 
used to advance the solution in successive time increments At. The process in- 
volved in using the finite-element method and the present timewise solution 
procedure follows (see the information flow chart of Fig. 4a) : 

Step 1 : Construct the mass matrix [ml for each finite element and then 

assemble these contributions according to Eq. 3.23f to form the mass 
matrix [M] for the complete assembled discretized structure. This 
[M] represents the "final" mass matrix if the structure has none of 
its generalized displacements constrained (that is, held equal to zero, 
for example) ; however, if such constraints exist, one forms a reduced 
or constrained mass matrix (and, in fact, a reduced set of the equations 
of motion) by deleting the rows and columns of [M] associated with those 
generalized displacements which are prescribed to be zero. 


Next, this constrained mass matrix is factorized to consist of a lower 

T 

triangular matrix [L] and an upper triangular matrix [L] according to 
the Choleski scheme: 


[L] [L]"" 


Since [M] does not change in value with time as the transient structural 

T 

response proceeds, one needs to determine [L] and [L] only once — these 
quantities need not be re-evaluated at each time step of the calculation. 

Step 2 : The prescribed externally-applied transient forces can be em- 

ployed to calculate the generalized applied forces {f} acting on each 

discrete element at each time instant t of interest. These, in turn, 

m 

can be assembled according to Eq. 3.16c to form the assembled applied- 
loads vector {f*} for the complete assembled discretized structure. 

Step 3 : Assuming that at zero time (t *> 0) , the generalized displace- 
ments {q*} = 0, the generalized velocities are nonzero {q*} = {a}, and 

o o 

that nonzero external forces {f*} q are present. In this case, Eq. 3.18 
becomes 
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[Mj \f\ = fF*j„ 

fL][L] T /Fi = IF\ 


from which one can calculate {q * } q by using the earlier-described Choleski 

scheme. Then from Eq. 3.29 one obtains 

. .2 . ^ 


toV tJ r 


(3.37) 


»«?*), -nt 

m- fF = prescribed initial generalized 


(3.38) 


velocities 


Also, 


ffri - sri -wi, 


(3.39) 


For this case, however, it has been assumed that {q*} = (o). Thus the dis- 

o 

placement configuration {q*}, at time t, = t + At is known. 

1 1 o 

Step 4 : Knowing the generalized nodal displacement increments 

{Aq*}^ = {q*> 1 - {q*} o ana the generalizea nodal displacements { q* } ^ 
at time t^, one knows also the unstarred individual element quantities 
{Aq}^ and {q}.^ via Eq. 3.15. Hence, one can calculate the strain in- 
crement (Ay. . ), developed from time t to t, at every Gaussian station 
i) 1 o 1 J 

(or point) required over and depthwise through each finite element 
from Eq. 3.10: 

( *y, -p, = )„ 

= ^f\i D 4 LD i J l3 - 40> 
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With a knowledge of (a) the stresses at t 


t^ - At, and (b) the strain 


increment (Ay^ J one can determine the stress increments (As 1 " 1 )^ 


and 


the stresses (S J ) ^ at time t^ at each Gaussian station by using the 
pertinent elastic-plastic stress-strain relations, including the 
yield condition and flow rule (this matter is discussed in detail in 
Subsection 3.3.2) . 

Step 5; Next, one can calculate {p}^ and [h] ^ for each individual 

finite element by using Eqs. 3.14a and 3.14b, respectively. Assembly 
of this information according to Eqs. 3.16a and 3.16b, respectively, 
provides {p}^ and [H] . Since the prescribed generalized force vec- 
tor {F*}^ is available from known {f}^ information, the equation of 
motion, Eq. 3.18, at time instant t^ becomes: 


CMJ^VfFVfP), -fHJ, fy'i, 

In the interest of minimizing computer storage and the number of 
manipulations, one first forms for each individual element {b R }^ = 
({f} - {p} - [h] {q})^. Then one forms the right-hand side vector 
of Eq. 3.41 by 

L J J T { ki , k.i, , ku j , 

For clarity of discussion, however, the form of the equation repre- 
sented by Eq. 3.41 is used here. 

Step 6 : Since the right-hand side of Eq. 3.41 is now known, one can use 

the Choleski scheme to solve the following equation for the accelera- 
tion {q*}^: 


(3.41) 


(3.42) 


(3.42a) 


Step 7 : 


With {q*}^ now known, one can calculate the generalized displace- 


ment increment {Aq*}^ from Eq. 3.25a as 


(3.42b) 
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where 


i-rt - \r \ - » A 

(3.43a) 

= f?V ffl 

(3.43b) 

Thus, from Eq. 3.43a one has 


I f \ 2 = i A + 

(3.44) 


The process then proceeds cyclically from Step 4 onwards for as many time 
steps as desired. 

For conciseness, the selection of a pertinent time increment size At 
is discussed in Subsection 5.3. 

3. 3. 1.2 Conventional Formulation 

For convenience, the equations of motion which cure to be solved in the 
conventional formulation (Eq. 3.24) are repeated here: 

(M] f f j * ( K ] \f 1 = f F'j f F‘ p * I F + ) * I F ^ (3.24! 

The solution process for Eq. 3.24 (see Fig. 4b) is very similar to that just 
described for the improved formulation except for some modifications in 
Steps 1, 4, 5, and 6. To avoid needless repetition, only these modifications 
are described here. 

Modifications to Step 1 : In addition to forming and factoring IM] , one 

must form the stiffness matrix [k] for each element and then assemble 
this information according to Eq. 3.23a to form [K] for the complete 
assembled discretized structure. Note that [M] and [K] need to be 
formed only once — they do not change as the transient structural 
response proceeds in time. As in the case of [M] , if displacement 
constraints are present in the problem being analyzed, one must 
form a "constrained" or "reduced" stiffness .matrix for the entire 
assembled discretized structure. 
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Modifications to Step 4 ; In addition to requiring the determination and 
storage of stress increments and stresses at each Gaussian station, 
one is required now to determine and store also the plastic strain 
increments and plastic strains at each Gaussian station at each time 
instant. 

Modifications to Step 5 ; Because of the presence of large deflections 

and elastic-plastic effects, one must calculate {f NI *} , {f 1 "}, , and 
NL q i P i 

{f^ for each finite element by performing numerically the volume 

integrals indicated by Eqs. 3.22, 3.20e, and 3.20g, respectively. 

The assembly of this information according to Eqs. 3.23b through 

3.23d is then accomplished to form and {f*^}^. 

Here again the actual operations are done more compactly than this 

description implies . 

Modifications to Step 6 : Now at time t^ all information needed to 

write Eq. 3.24 is available: 

[ uoif f\ -tKH (fF'J - f Fp * i f; l ] - 1 Fti ) , , 3 . 45) 

Note that the forming of the right-hand side (RHS) of Eq. 3.45 
requires the multiplication of [K] by {q*}^; then this “force 
vector" is added to the remaining terms of the RHS force vector. 

Thus one needs to have both [LJ and [K] stored for use. 

Otherwise the timewise solution process for the conventional formu- 
lation proceeds in the same fashion as that described for the im- 
proved formulation. 

3.3.2 Evaluation of Stress Increments, Stresses, Plastic 
Strain Increments, and Plastic Strains* 

In the present subsection, the calculation procedure for the determina- 
tion of stress increments and stresses at any station (such as Gaussian, for 
example) in each element is described. Because the "mechanical sublayer model" 

+ 

This procedure applies for both the conventional and the improved formulation. 
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is adopted in the present analysis, the only constitutive relation considered 
is that for a homogeneous , initially isotropic, elastic, perfectly-plastic 
strain-rate dependent solid; strain hardening is automatically accommodated by 
the mechanical sublayer model. 


A convenient way to compute the stress component increments and stress 

components at time t = (m+1) At, as discussed in Refs. 46 and 48 will be 
m+x 

employed? it is assumed that all stresses and strains are known at time t . 

m 

One begins by assuming that the strain increment (Ay . , ) , , , from time t to 

13 m+l m 

time t ^ as calculated by Eq. 3.40, is entirely elastic, and a trial (over- 
script T) value of stress increment is calculated from the relation* : 


( aS ■ ) 

K J 'm+l 


% [ L? 1 -zy (A if t J 


where 


x k 


ay; - 


Hence, the trial stresses at time sure given by 


(S> )„+» 


= cs;-i - (^5;- u, 


(3.46) 


(3.46a) 


(3.47) 


Then a check process is performed by substituting this trial value of the 
stress into the Mises-Hencky yield function, Eq. 2.79, to determine whether or 
not the trial stress state lies inside the yield surface; thus one may write 

; T - JtiVcM fCl ) 1 2 2 




(3.48) 


where is the appropriate known uniaxial yield stress. Also, it should be 

recalled that the Mises-Hencky yield function is expressed in Eq. 2.68 in terms 

of T* where T* * i/g/G S 1 .. 

3 3 3 

If 4> < 0, the trial stress state lies within the elastic domain 

m+l — 

bounded by the yield surface. Therefore, for this time step there has been no 

plastic flow and the incremental deformation can only be elastic. Hence, the 

actual, stress (S 1 ) . . is equal to the trial stress; thus 
3 m+x 

T . 

(5J ) w+/ = ( Sf (3.49) 

+ 

Such calculations are carried out for each layer of the mechanical sublayer model. 
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and , the plastic strain state is 


j, ■ 

( ) = ( -y J ) 

K * 1 '* 11 + 1 ' f m 


(3.50) 


However, if $ > 0, the trial stress state lies outside the yield surface 

m+i 

(i.e. , in the undefined region). Therefore, the trial assumption that the 
entire strain increment in an elastic-strain increment is not valid. Plastic 
flow has occurred within this time step and the actual stress state must lie 
on the yield surface according to the theory of perfect plasticity. Then, 
the calculation proceeds as follows. 

The total strain increment can be decomposed into elastic and plastic 
components 




<*y}L, = < 4Y p-*' 


(3.51) 


The stress increment is related to the elastic component of the strain 
increment by the relation 


( 4 s ' + ; & A ^i. s i\,i 


and the actual stress is 




(3.52) 


(3.53) 


Since the material is assumed to be incompressible with regard to plasticity, it 

■b Jq ^ ^ 

follows that Ay, «* 0 (orAy . = 0: spherical component of the incremental plastic 
* 3 

strain). Then from Eg. 3.51 one has 


and 




(3.54) 


(3.55) 


Substituting Egs. 3.54 and 3.55 into Eg. 3.52, results in 
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* y i * r £ r f - 


(3.56) 


where Ay is given, in general, by the plastic flow rule, Eq. 2.82, as 

Vf ■ & ■ 

A ~^ 7 “ A ( 3 . 57 ) 

where 1 is a real nonnegative constant which will be calculated later. 


The stress S_. determines the direction (or relative proportions) of the 
plastic strain increment. These directions will vary, over the time interval 

°i D i 

At from (S.) to (S.) as a result of continually straining. However, for 
3 m 3 m+1 &T ^ 

computational convenience, the deviatoric component of the trial stress (S . ) , 

3 »+l 

which lies between (S.) and (S.) . , will be used to approximate the correct 
3 m j m +1 

direction. Thus, Eg. 3.56 becomes 


( A Sfh.i i+y 




(3.58) 


Further, combining with Eg. 3.46, one has 


^m+i ~ ^ A Sj K+t ^ S j )„*, ^ 


where 


^*,*1 ~ \+ i > ' 


Then, the actual stress at time t , is 

m+1 




(3.59) 


(3.60) 


The plastic strain at time t , is given by 

m +1 


( ~ ( ( m + ( i*' ? ' h,+ 1 (3.61) 

The quantities X^ + ^ and X! ^ in Eqs. 3.60 and 3.61 can be determined 
from the fact that ( S j) m+ j_ must satisfy the yield condition: 


Alternate approximations could be used (see Refs. 133 and 134, for example). 
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one obtains the 


Substituting Eq. 3.60 into Eq. 3.62 and solving for X ,, 

SH*J. 

physically valid values 

c 

7\ 


where 


w+l B +] B*~A C 

J?T • J>T • / DT . 2 

A - ( 5 ) ( 5 \ ~ kJm.t 


L T 


b.(s;u. 


vtiri • — k J m+- 1 


'" T ‘- csi). 


G . 


2 6 


f - ( siijsiLrs^tL - j ? r / 


(3.63) 


(3.63a) 


The preceding discussion has pertained to the use of elastic, perfectly- 
plastic rate-independent material whose yield stress is 2 O q . However, if 
the yield stress is rate dependent, the same procedure applies except that the 
yield stress in Eq. 3.48 is the strain-rate dependent yield stress which is 


given approximately by Eq. 2.87 as 

<K = <£ f i + \jr\* J 


(2.87) 


where O is the static uniaxial yield stress, D and p are material constants, 
o 

and e is the uniaxial strain rate. For the three-dimensional problem, it is 
assumed that e of Eq. 2.87 may be replaced by the second invariant of the devi- 
atoric strain-rate tensor. Thus, this equivalent strain rate is given by 


e = J f ( fi - j < ~nV) 


(3.64) 


* i * j * i i 

where the strain-rate components Yj (and/or yi) are given by y_. = (AyJ/(At). 
Other alternatives for t have been proposed (see Refs. 46, 133, 134). 

It should be noted that in the solution procedure for large-deflection 
elastic-plastic transient responses, the tracing of the st--os history is 



required for both the conventional formulation, Eq. 3.1, and the improved forma- 
tion, Eq. 3.2, because of the nature of the elastic-plastic theory used. How- 
ever, in addition to the stress history, the tracing of the plastic strain 
history is also required for the conventional formulation but not for the im- 
proved formulation. This may result in the saving of computer storage, if the 
improved formulation is used. Further discussion of this will be given in the 
next subsection. 

3.3,3 Comparison of Storage and Computing Requirements for 
the Improved versus the Conventional Formulation 

From the solution process discussed in Subsections 3.3.1 and 3.3.2, it 
is clear that the storage/computing required are less for the improved formula- 
tion (1-method, for short) than for the conventional formulation (C-method) . 

By comparing the conventional formulation, Eq. 3.1, with the improved formula- 
tion, Eq. 3.2, it is seen that the storage of both the assembled mass matrix 
and the assembled stiffness matrix are required by the C-method but only the 
assembled mass matrix is required by the I -method. Also at each time step, 
the matrix multiplication [K]{q*} is needed for the C-method but not for the 
I -method. 

Further, as was mentioned before, because the structure undergoes large- 
deflection elastic-plastic behavior, Gaussian integration is employed to evalu- 
ate {p} and [h] in the I-method; hence, the storage of the stress history at 
every Gaussian station in each discrete element is required by the I-method. 

On the other hand, if the C-method is used and Gaussian integration is also 

employed to evaluate {*“■}. {f 1 } and {f NL }, then, in addition to the storage 

*1 P P 

of the stress history, the storage of the plastic strain history at every 
Gaussian station in each discrete element is also required. As for the com- 
puter operations, at each time step, three matrices {f 1 ^}, and {f L }, and {f NL } 

*1 P P 

need to be evaluated for each discrete element, if the C-method is used. But 
only two matrices: {p} and [h] (or {p} and [h]{q}) need to be evaluated for 
each discrete element, if the I-method is employed. 

Based on the above storage and operation considerations, it may be con- 
cluded that the I-method is more efficient and simpler than the C-method. This 
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conclusion is substantiated by the computing experience discussed in Section 5 
for various illustrative examples. 

3.3.4 Selection of a Temporal Finite-Difference Operator 

For the timewise numerical solution of undamped linear dynamic structural 
problems, many finite-difference operators have been explored to assess their 
attributes and shortcomings. Some schemes are stable no matter how large the 
time increment At is chosen to be — and hence are termed "unconditionally 
stable"; others are unstable for At larger than some critical value — and 
thus are termed "conditionally stable". Some introduce (unintentionally) arti- 
ficial or false damping whereas others do not exhibit this undesirable feature. 
All of these methods, however, usually* produce a phase-shift error in the pre- 
dicted response, depending upon the size of the finite At used — some schemes 
exhibit more phase-shift error than others for a given At. A concise tabula- 
tion of some features of the more commonly-used varieties of this method is 
given on the next page (Refs. 132, 135-142) together with some examples of 
users of each method for linear and/or nonlinear transient response predictions. 

The criteria for stability of each of these common methods have been 
established for linear transient response problems (Refs. 143-151) . These 
studies have derived the At conditions under which exponential round-off error 
growth will result. For smaller At values, this type of error growth will not 
appear, and a "stable" calculation is said to result. O'Brien et al. [143] , 
Leech et al. [146, 147, 148], Johnson [149], Nickell [150], and Krieg [151] 
have illustrated this type of analysis and behavior. Nickell's study [150] is 
especially extensive, treating the 3-point central-difference method, the 
Newmark 6-method, the Wilson averaging method, and the Gurtin averaging method. 

It should be noted (Ref. 48, for example) that one can readily construct 
m-point forward-difference, central-difference , or backward-difference opera- 
tors by Taylor series representation of the accelerations x and/or velocities x 
* 

An exception has been noted in Ref. 146 wherein the 3-point central-difference 
formula was used to solve the one-dimensional wave equation. When At was chosen 
such that (At) /(Ax) = 1, a solution which was exact in both amplitude and phase 
was obtained. 
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in terms of displacement information at m instants in time; the truncation 
error of each approximation thus selected may be readily identified, and de- 
pends upon the number (m, such as 1, 2, 3, 4, etc.) of time instants used. 

Further by using the methods of Refs. 143, 145, 148, and/or 150, it can readily 
be shown that: (1) all forward-difference operators are unconditionally un- 
stable , (2) all central difference operators are conditionally stable (a criti- 
cal At exists beyond which error blowup will occur) , and (3) all backward- 
difference operators are unconditionally stable (i.e., stable for all At). The 
Houbolt method is a four-point implicit backward-difference method (that is, at 

time n, x and x are expressed in terms of x , x . , x and x ,) ; this 

n n n n-1 n-2 n-3 

method accordingly is unconditionally stable [149] . Note that all of these 
implicit methods cited in the tabulation on page 70 are unconditionally stable 
except for type 12(b) — a version of the Newmark 2-method. Methods 12 
through 15 were not constructed from the above-described Taylor's stress approxi- 
mation — somewhat different intuitive and/or rational procedures were used. 

Note that all of the implicit methods except the 2 = 1/4 version of 
Newmark' s 2-method introduce false damping. The latter method and the 3-point 
central-difference method noted in the tabular summary do not introduce false 
damping. In Newmark 's 2-method, for example, the amount of false damping de- 
pends upon the value of 2 used; Newmark suggests (a) choosing 2 - 1/12 if one 
seeks high prediction accuracy for an extended period of response for a struc- 
ture with small actual damping or (b) choosing 2 = 1/6 if one is interested in 
only a few cycles of response — the implication being that the error introduced 
by false damping would be acceptably small for many engineering purposes. 

While round-off error instability is avoided by all of the unconditionally 
stable methods (permitting one to use as large a At as one wishes) , the forcing 
function in a given problem may have severe spatial and temporal variations such 
that one must use a fairly small At in order to follow and identify the severe 
peaks, etc. in the structural response. Perhaps a At of some chosen fraction of the 
period of the highest significantly-excited mode should be used — provided that 
one can make a rational estimate of this situation. In such cases, the feature 
of unconditional stability may not be as much of an advantage over a conditionally 
stable method as one might think at first sight. However, for transient loadings 
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which excite only the lowest few modes of the structure, the "larger At" per- 
mitted by the unconditionally stable methods compared with the "stringently 
small At" permitted by the 3-point central-difference conditionally stable 
method can be clearly advantageous. 

Although one can construct finite-difference operators of the imp licit 
type or of the explicit type having truncation errors as small as one wishes 
by using information at time stations (n,n-l, n-2, n-3, ... ) or (n+1, n, n-1, 
n-2, n-3, ... }, respectively, it is evident that one pays a price in the 
necessity of storing this information in order to march the solution ahead in 
time. Further, for large-deflection problems involving elastic-plastic be- 
havior, the use of an explicit operator circumvents the iterative type of 
calculation (or extrapolation) for the equivalent generalized loads required when 
an implicit operator is used. These considerations indicate that one should 
choose an explicit operator whose accuracy vs. storage tradeoff is most bene- 
ficial. In view of its simplicity, accuracy, lack of false damping, and mini- 
mal storage required, the 3-point explicit central-difference operator has 
been chosen for principal use in this study; studies to define an "optimum” 
operator of this type have not been carried out. 

Although the criteria for stability of each of these common methods have 
been established for linear transient response problems (Refs. 143-151) , no 
similar assessment is known to have been made when these methods are applied 
to nonlinear structural response problems involving large deflections and in- 
elastic material behavior. Various of these methods, however, have been 
applied to such problems — with At values chosen in conformity with the estab- 
lished stability and/or convergence criteria for their use on linear problems, 
or by numerical experimentation. 

It has been demonstrated in the present study that the (a) Houbolt and 
(b) Newmark ((3 = 1/4) method both of which are unconditionally stable for linear 
structural response problems, now both become conditionally stable for large- 
deflection nonlinear responses whether the material behavior is linear elastic 
or elastic-plastic (see Subsection 5.3.2); also the 3-point central -difference 
method remains conditionally stable but the stability criterion becomes more 
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severe (smaller At is required) them for linear problems . These results con- 
firm similar findings by Stricklin [75] and McLaughlin [156] . Further dis- 
cussion of this matter is given in Subsection 5.3.2. 
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SECTION 4 


FORMULATION FOR A CURVED BEAM 


4 . 1 Introduction 

In Section 2, the equations which govern the large-deflection, elastic- 
plastic dynamic responses of a general 3-dimensional continuum are described. 
Section 3 presents the overall method of analysis. Based on the Principle of 
Virtual Work and D'Alembert's Principle, the spatial finite-element approxima- 
tion has been used to derive the equations of dynamic equilibrium. Then, a 
direct numerical integration scheme with an appropriate timewise finite-dif- 
ference approximation is used to solve the resulting equations of motion. 

In the present section, the application of this approach is demonstrated 
in detail for curved beamlike structures which undergo planar (two-dimensional) 
deformation with or without the inclusion of transverse shear deformation ef- 
fects. In the structural finite-element context, such configurations are 
termed "one dimensional". 

An arbitrarily-curved beam element is described here. Its specializa- 
tion to represent simple circular ring and straight beam structural elements 
is given in Appendix B. 

The geometry of a curved-beam element is described in Subsection 4.2. 

The formulation for a Bernoulli-Euler-type curved beam element is discussed in 
Subsection 4.3, while Subsection 4.4 is devoted to a corresponding development 
for a Timoshenko- type (shear deformable) curved beam element. Both small- and 
large-deflection behavior are included. 


Description for a Curved Beam Element 


The geometry and nomenclature of a typical undeformed curved beam element 
(Refs. 157 and 158) are shown in Fig. 5. The parametric equation of the beam's 
centroidal axis on the planar surface can be expressed as 


n = n <? ,= - z»t> k 
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where r) is the length coordinate measured from node i along the centroidal axis 
(meridian) and XYZ represent global reference Cartesian coordinates. 


The unit tangent vector, a, to the centroidal axis, and the unit normal 
vector, n, cure defined as 

_ 2 _ 

— ]_ d a. _ _ d n 

dr\ . n = > dr( pi *} 1 


*<>?> - dn 


(4.2) 


where p is the magnitude of the curvature vector da/dr), and the reciprocal of 
p is known as the radius of curvature, R, and is taken positive when the center of 
of curvature lies in the negative direction of n. 


It is assumed that the slope, of the centroidal axis, which is the 
angle between the unit tangent vector and the y-axis of the local reference 
Cartesian coordinate system (x,y,z) may be approximated with sufficient accuracy 
by a second-order polynomial in n as follows: 

4-pp = b. + b, >z + (4 . 3) 

The constants b Q , b^, and b 2 can be determined from the known initial geometry 
of the curved-beam element by requiring (1) the slopes of the idealized ap- 
proximated beam element and the actual beam element to have the same slopes at 
nodes i and i+1 and (2) the ends to lie on the y-axis (i.e., z = 0 at both 
ends) . Thus 

<4 < 7 = 7* = 0 ) = (4.4) 


<M ? = ) = 4> xW 


(4.5) 


ft" &N<t> d*l = 0 

'o 

If <j> is small, Eq. 4.6 may be approximated by 

4 d 7 = o 

Using Eqs. 4.3, 4.4, 4.5, and 4.7, one obtains 


(4.6) 


(4.7) 
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b, 


b i 

b z 


= 

-Z (4^, -W/ ) 

7/* i 

_ 3 ( i, 4 <k ) 


(4.8) 


Accordingly , 
pressed as 


the radius of curvature, R, of the centroidal axis may be ex- 



Consider the beam subjected to 2-dimensional deformation. A generic 
point p in the beam element is displaced to a new position P. Its new position 
vector, R, is given by 



r + V 


(4.9) 


where r is the position vector to point p, v is the displacement vector which 

is a function of (n , C) and £ is the distance of point p from the centroidal 

axis along the unit outward normal, n, direction. The displacement vector v 

n, 'v, 

may be written m terms of its components denoted by v and w in the direction 
of a and n, respectively. Thus 

7 = 9" ( ? . s ) a + yv (>?, £ ) n 

(4.10) 


and the displacement v(n,C) and w(n,C) may generally be expanded in power 
series of £ by 


"v<7- ?)= + ’ 




w 


(4.10a) 


4.3 The Bernoulli-Euler-Type Curved Beam Element 
4.3.1 Displacement Field 

Let the Bernoulli-Euler hypothesis (Refs. 11 and 159, for example) that 
the beam cross section which is perpendicular to the centroidal axis prior to 
the deformation remains plane and perpendicular to the deformed centroidal locus 
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after deformation, and also that it suffers no strain in the direction normal 

*\r 

to the centroidal axis be employed. Hence, let the displacement field v,w of 
the beam, Eq. 4.10a, be approximated by the middle plane displacement v and w, 
and the rotation ^ as follows: 

7 (7 . s ) = vc%)- i; ^(?) 


where 


=ff ~£ 


(4.11) 


(4.11a) 


The selection of a suitable interpolation function to represent each of 
those displacements throughout each element is one of the principal concerns 
in the construction of finite-element assemblage of the whole structure. It 
has been shown in Ref. 160 that the inclusion of rigid-body-displacement modes 
in the assumed displacement field of a cylindrical-shell element leads to a 
better coarse mesh solution than if the rigid-body modes are excluded from the 
assumed-displacement field for the linear-elastic static cylindrical shell prob- 
lem. Also, it has been concluded in Ref. 161 that the use of a cubic polynomial 
to express both the axial displacement v and the normal displacement w in the 
circular ring element exhibits a marked improvement over the use of a linear 
expression for v and a cubic expression for w; also, the former converges very 
rapidly to the exact linear elastic static solution. Based on these considera- 
tions, two sets of assumed displacement functions for the present Bernoulli-Euler- 
type curved beam element will be formulated in the present analysis: (1) both v 
and w will be represented by cubic polynomials in T) with rigid-body modes in- 
cluded (the finite-element formulation from this expression will be denoted as a 
CC, or cubic cubic, element) and (2) a linear expression in n for v and a cubic 
expression in n for w, also with the rigid body modes taken into account (this is 
termed the LC, or linear cubic element) . 


Assuming that the element is subjected to small amplitude rigid-body 

translations V and V , and rotation (2 with respect to the local reference 
y a x 

Cartesian coordinate system (x,y,z), the rigid-body displacement expressed along 
the curvilinear directions n, a of any point p(y,z) is given by 
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! V l 

u 


co5 4> 4 1 -zcos<$ +ysiH<p "j 

“SiN 4 s cos^ zsiN4>+y 


^y 

v * 

n t 


(4.12) 


r/jid body 

To account for the strain- inducing inodes and the rigid-body modes, the 
assumed displacement field for the CC element takes the form (Ref. 160) : 


C0S4> SlN<p -Zcosdf+ySMfy *1 0 0 Y] 1 >p ? j j # 

-SIN CDS 4 ZSItity+j C°S 4 0 Y} 1 0 0 


I 

h 


Pe 


(4.13) 


where 3 , 3 2 , 


, & are parameters which will shortly be expressed in terms 

O 


of the eight selected generalized displacements of the element. In more com- 
pact matrix form, Eq. 4.13 becomes 


(uj = [ U <?>) I 

1 - " 8, I 


(4.13a) 


2*1 z * 8 

The generalized displacements, termed {q}, are chosen to characterize the de- 
formation state of the element, and are selected such that there are four de- 
grees of freedom v, w, y, and v, (= 3v/3ri) at each node of the element: 


1 


= L 


vV,- 


% 


T 


v i+i 


t, 


i+i 


J 


(4.14) 


Corresponding to the assumed displacement field, Eq. 4.13, one finds 
_ JW v' _ „ „ , n, _ vi 2 . w3. 


15) 


0 I -?/R <?? 5? 

_ ay 

— -L~w ~fT Ft 1 y y c '( 3 7JTp; 

(4.16) 


y 


= L S “^ 2i-z™*+ysi£*l , 0 o 3 >7 

l 


'7 ^7 L R R <? T 

= L^k, ? j /p] 

The generalized nodal displacements, {q}, and the parameters, (3), of the 
assumed displacement field are related by a transformation matrix [A] which 
may be obtained by substituting the coordinates of nodes i and i+1 into Eqs 
4.13, 4.15, and 4.16. Thus 


m = (4) foj 

8*1 8 <8 8*1 

Because [A] is a square nonsingular matrix, one may write: 


(4.17) 
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(4.18) 


IJ 3 ) - ( 

S*I 8*0 8*1 

Substituting Eq. 4.18 into Eq. 4.13, one obtains 

U j = [Lh^JM) \f\ (4.19) 

Z.x I 2 k3 8x8 8*1 


Turning now to the LC formulation, the assumed displacement field takes 
the form 

I ft 


j v l _ r cos4> sin 4> -zcos4>+ysiN<t> y 


IwJ 


1-VSIN4* costp z sin< 1> +y cos<f> o +i l 


k 


(4.20) 


and the generalized nodal displacements now chosen to characterize the deforma- 
tion state of this LC element are defined to be v, w, and ip at each end of the 
element. Thus 


{ =L W, Vtvi v^v, (4 ‘ 21) 

6*1 6*6 6*1 

It should be noted that by the nature of the assumed displacement 
finite-element variational principle used, the internal equilibrium equation as 
well as the force boundary conditions are generally not satisfied everywhere 
exactly by its solution, although the displacements obtained by this method are 
usually a very good approximation. The effect of using more displacement modes 
is to improve the satisfaction of the equilibrium condition in the interior of 
each individual element and hence also the accuracy of the approximate solution. 
However, the compatibility of the additional displacement mode, such as v,^ 
in Eq. 4.14, with neighboring elements is not necessarily required from the 
point of view of defining and evaluating the variational argument in Eq. 3.3 
because the strain depends only on the first derivative of v with respect to n, 
as will be seen in the next subsection . For a static problem, the two gen- 
eralized coordinates v and v,^ + ^ of each element may be condensed out 
through the use of the static condensation process (Ref. 7) , but a rational 
condensation process has not yet been devised for the corresponding dynamic 
problem. In the present CC-type element, the v will be treated as an 
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independent generalized coordinate at each node of the finite-element assemblage 
and the compatibility of v, between neighboring elements is retained. 


In the following, the CC-type assumed displacement field is employed in 
conjunction with the strain-displacement relations to obtain the equations of 
dynamic equilibrium from the pertinent variation statement, Eq. 3.3. The 
finite-element properties and the governing equations for the LC-type assumed 
displacement field may be derived in a similar manner, except that the corres- 
ponding fUJ _ , [A] matrix should be used. The various matrices which are 
zXo bxb 

used symbolically in this section of the text are documented in Appendix B. 


4.3.2 Strain-Displacement Relations 

Under the Bernoulli-Euler hypothesis , the only nonvanishing stress com- 
ponent and corresponding strain component which need be introduced into the ap- 
propriate beam theory are the axial stress 0 and axial strain, e. Also, let the 
ratio of thickness to radius of curvature be negligible in comparison with unity 
The appropriate nonlinear strain-displacement relation may be expressed as 
(Refs. 9 and 159) : 


where 


£ = E* 


C K 


= 


*7 


vV 

R 


K = 


R > 

?W - _ <L f 3 w _ X 


_L i ?w _ x. ' z 

z 1 


or in matrix form: 


E„ = L 


/ 

I 


i>7 ( 3 7 


3*7 R-l 


( V l 




n 


(4.22) 


(4.23) 


w 


} + IL V ^ £ 


2- , J v L 
L R ?7 J ) w J 


= l8,j {«} ^ y L<i_l { 6 2 j l S* j , 5 iij 


K -Lifi?*} ly) pfu j * j = l j 


(4.23a) 


Combining Eq. 4.19 with Eq. 4.23a, one obtains 
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(4.24) 


e. = - i LfJ {D t \mi \%\ 

k = lBjjI?! 


where 

lD.j 


l 6,j t/n 


In the process of solution, it is necessary to evaluate the strain incre- 
ment, Ae^, from time t^ to time t^. Using Eqs. 3.40 and 4.24, the strain 
increment is related to both the displacement and the displacement increment by 


- * t ai + £ A K a 


(4.25) 


where 

^e a< = LC, j \*f\. T l "Ki d 4 lD ^ j ^ ^ 


(4.25a) 


^ k; = l j { A ^ | A - 

(4.25b) 


In Eq. 4.25a, {AqK is the generalized nodal displacement increment from time 
t^_^ to t^, which is computed directly from the equation of dynamic equilibrium 
of the system, {qK is the generalized nodal displacement at time t^ and is 
{q}^ = {q}^_^ + {AqK. The last term in Eq. 4.25a is of higher order compared 
with the other two terms, if the time increment step. At, is small. However, 
this term can become significant for the case of a sufficiently large time incre- 
ment step. Its effect will be discussed later in Section 5 in the context of 
numerical examples. 

It should be noted that the only nonlinear term retained in the strain- 
displacement relations, Eq. 4.23, is due to the rotation of the centroidal axis. 
This expression for the strain is suitable for the cases where the deflection is 
large compared with the thickness of the beam, but it is still small compared 
with the spanwise (longitudinal) dimension of the beam. Otherwise, the follow- 
ing more accurate displacement and strain-displacement relations should be 
used (Ref. 9) : 
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(4.26) 


and 


Vl C ) = V(>?) - 
w (?, r ; ) = vv ( 7 ; + £ -x ) 


where 


V = 


31V 

*7 


x 

R 

w 


c — ( ^ jd 4 . — ) 4 . 

£ ° - ( J7 R ; 


-v - IX 4 *L 

^ - dr/ R 


1 ( 2Ji 

2 1 17 


v / +. 1 ( 2X + w 1 1 
' 2 ( a? R ' 


R 


*. - H + r-Vt +2. 
K - ( 1 R ) ( ?7 - 


3»v 


R } ^ ( ^7 ~ R U 2? 


?Z +JP 


(4.27) 


) 


(4.27a) 


However, in the numerical examples carried out in the present analysis, only 
Eqs. 4.11 and 4.23 are used. 

Also, it should be noted that the rigid-body displacements given by 
Eq. 4.12 are only approximate because of the assumption that the amplitude of 
the rigid-body rotation ( 2 ^ is small; this displacement field yields zero 
strain when applied in the strain-displacement relation given by Eq. 4.23 (or 
Eq. 4.27) only for small deflections. If the element is subjected to large 
amplitude rigid-body translations V^, and V^, and rotation with respect to 
the (x,y,z) coordinate system, then the correct rigid-body displacements ex- 
pressed along the curvilinear direction a,n of any point p(y,z) would be 

j V t _ I COS 4> S/N^l j V y l + 

' wj (~SIn4> Cos4>J { J 


} A C.OS <p ■+ BSIN <P 
-A sitf 4> + B cos dp 


(4.28) 


where 


A = y cos - z sin fl x ~ y 
B = y sinJ7 x cos SI x -z 


(4.28a) 


These large amplitude rigid-body displacements expressed by Eq. 4.28, will re- 
sult in the prediction of zero strain when the more accurate strain-displacement 
relation Eq. 4.27 is used, but not when the approximate strain-displacement rela- 
tion represented by Eq. 4.23 is applied. However, in the present analysis, the 
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"simplified" strain-displacement relation Eq. 4.23 together with the small ampli- 
tude rigid-body displacement modes represented by Eq. 4.12 are employed prin- 
cipally because of the unwieldy nature of the expressions resulting from using 
Eq. 4.27 together with an assumed displacement field which includes both the 
proper deformations and the exact rigid-body displacements given by Eq. 4.28. 


4.3.3 The Discrete-Element Property Matrices 

Under the Bernoulli-Euler assumption, the consistent mass matrix of the 
discrete element including both rotatory and translational inertia effects may 
be obtained from the expression for the kinetic energy, K^, as follows: 

K 


or 


where 


1C = 


= i ((( Po dv n 




= 7 ff( Pol ( v-^) 2 - 

. 2 

w J 

dV n 

(4.29) 

= 7 L * w ^ j f B ] 

1*1 
l -y ; 

dy l 

(4.29a) 

r K o 

0 

1 



l B ] = 


0 

0 


P'b h 0 

o " ?±kK 
12 


(4.29b) 


and, b is the width, h is the thickness of the beam, and p is the mass density 

o 

per unit volume of the undeformed body. 

With the assumption that the velocity field is of a form which is con- 
sistent with the displacement function, Eqs. 4.13 and 4.15, one has 



= [ N<7'] I f\ 

3*8 8 ’ i I 


(4.30) 
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wnere , n 

l\(>?) -i 

r M c^j] = - - - - 2 - 8 [A] 

1 J J G^vi 8<i (4 - 30a) 

1 1 D 

is consistent with Egs. 4.13, 4.15, and 4.18. 

Substituting Eg. 4.30 into Eg. 4.29a, one obtains: 

K E = yLjj|^'[N(?)J T [B] IN <?>]«*? If] (4.3i) 

or 

K e = L r m 3 ^ ^ J (4.31a) 

where the consistent mass matrix [m] of the element is 

P") =f V '"[N(V 1 ] , (6)[N<'?>] rf? 

■>v 

The element's consistent mass matrix defined by Eg. 4.31b is symmetric 
and positive definite, and is consistent (in the variational sense) with the 
assumed displacement field. The resulting eguations of motion obtained by 
using the consistent mass matrix will be a system of simultaneous (coupled) 
second-order ordinary differential eguations. In order to take advantage of 
potential storage and computing efficiencies, another mass matrix called the 
lumped (diagonal) mass matrix (Refs. 68, 73, 153, 162, 163) is often used in 
dynamic analyses; the resulting eguations of motion will be a system of de- 
coupled second-order ordinary differential eguations. The lumped element mass 
matrix of the present beam element can be written as 

r*S ~\ 


CH = p B b h 

lumped 


'5s 3 

"Ss 3 


(4.32) 
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where s = = element size, a and 6 are constants, and taken to be a = 1/2, 

6 = 1/24 (Ref. 162) ; other values also have been chosen for the constants a and 6, 
Ref. 73. However, further studies are needed to develop appropriate lumping for each 
of the various types of selected criteria such as (a) by frequency matching in which 
the lumped mass properties are chosen such that the lumped-mass system and the con- 
sistent-mass system have the same highest natural frequency (Ref. 73) ; (b) by stati- 
cally-equivalent considerations wherein the lumped mass properties are chosen to be 
statically equivalent to the actual mass distribution (Refs. 68, 153, 162, and 163) , 
etc. Also, some very useful information for this type of analysis is described in 
Ref. 190 in which the rates of convergence of the mode shapes and frequencies by 
the finite-element method using consistent- and lumped-mass formulations are 
established. 


The equivalent generalized nodal forces which correspond to or represent 
the externally-applied loading can be obtained by placing the assumed displace- 
ment field into the expression for the variation of the work of the externally- 
applied loading: 

§W = 

where 

F(t) = F v (t)a + F^(t) n is the applied time varying force per 

unit length 

M(t) = M(t)i is the applied time-varying moment per 

unit length. 

Substituting the assumed displacement function, Eqs . 4.13, 4.15, and 
4.18 into Eq. 4.33 


f 7,r '( p tt) S’ v + F„(t> Sw + M it) Sv ) d*i 
he v 


(4.33) 


— /%+* T 1 Fy 1 

Sw = Jf-i r n ( *?>] f F> j dy 


where 




if s 


1 T I 


m= w] p- 

J 7, Im 


(d = generalized nodal force 
matrix for the element 


(4.34) 


(4.34a) 


The equivalent nodal force which corresponds to the internal axial stress , 
0, also can be obtained by placing the assumed displacement field into the ex- 
pression of the variation of the work of the axial stress : 

85 



(4.35) 


£U=fff o-Zed v = fff cr(Ze B +$Sk ) a V n 

Vn V n 

Substituting Eg. 4.24 into Eq. 4.35 and introducing the stress resultants for 
the beam cross section 

L=/T<T<M , M =/(<’'? J A ,4.36, 

A A 

where the integrations being taken over the cross section. A, of the beam ele- 
ment, L, is the internal force, and M is the internal bending moment of the 
cross section, results in 


JU - L JM/^( {«jL*{qjMy? fhjuuL'? m 

= L Zp { \f\ + [M i fr\ ) 

(4.37) 


where y. 

If 1 = pf jhjL -fpM )i 7 

th) - /£ 


(4.37a) 


The integrations along the centroidal axis length of the beam element which 
appear in {p} and [h] of Eq. 4.37a may be performed numerically, for example, 
by using the Gaussian quadrature scheme (see Ref. 131) . The axial force L 
and moment M at those spanwise integration stations will be described and 
evaluated in the next subsection. Note that {p} and [h] are quantities perti- 
nent to the improved formulation . 


In the conventional formulation , the variation of the work of the axial 
stress, 6u, is expressed in terms of displacements, and the plasticity effects 
are taken into account through the use of "effective plastic loading”. This 
formulation will be described in the following: 


By substituting the relation 


86 



<r = E C£ - e. p ) = E (£ 0 c k - ^ ) 


(4.38) 


into Eq. 4.35 one has 

${J = fff E ( e 0 + £ k -e* ) c % + g $ * ) <J V n (4.39) 

V, 

Employing the strain-displacement relation, Eq. 4.24, Eq. 4.39 becomes 

SU = lS gj d'"< jD.jeU lD,j* j D s j !gj 

A U J\ u 

+ Ebh (T L FiM LD ' J ty^ i 1 v \ ir L 

+ /£' EU ( +f L ri c <i LD 'J^i ) f D 4 LD ' J ^ 1 ?!) 



(4.40) 


(4.40a) 


i f, L ! = /( j?" ( E £ : V' s + ?E:& f jh, J ) dV n 

A U 
A « 
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i f j L j = -</£' e kh < i ^p N ip >!^p? 

^ /£' E b b (L^jifi +r L f J Jjhj^-^j) <«■««> 

In Egs. 4.40b and 4.40c, is the total plastic strain at the end of 
the ith time step. Thus 

h i-l P „ c P 

£ f =H (4.41) 

A mVT m 

and the integrations along the length of the beam element are also performed 
numerically. The plastic strain increment, A E.^, and the integration of the 
total plastic strain e"T over the cross section, A, of the beam at those span- 
wise integration stations will be described next. 

4.3.4 Stress-Strain Relations 

Because of nonlinear material behavior, although the strain variation 
through the beam thickness, by the Bernoulli-Euler hypothesis, is linear, the 
variation of stress across the thickness may be nonlinear. For computational 
convenience, the stresses are evaluated at selected Gaussian points across the 
thickness, and the corresponding weighting factors are used in evaluating the 
pertinent integrals by Gaussian quadrature. The strain-hardening behavior of 
the material may be accounted for by using the mechanical sublayer model in 
which the material at each Gaussian station is treated as consisting of equally- 
strained sublayers of elastic, perfectly-plastic material, with each sublayer 
having the same elastic modulus but an appropriately different yield stress, 
as described in Appendix A. 

It should be noted that within the framework of the Bernoulli-Euler 
beam theory, although the transverse shear strain y is zero, the transverse 
shear stress T, is nonzero. With the presence of both axial stress a, and 
transverse shear stress T , the Ilises-Hencky yield function may be written as 

§ = - <rj (4 .„, 

where is the yield stress of the idealized elastic, perfectly-plastic kth 
sublayer in the uniaxial-tension engineering stress-strain diagram (see 
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Appendix A}. The process presented in Subsection 3.3.2 of calculating stresses 

and plastic strains will not function properly, if the yield function (Eq. 4.42) 

is used. (Because in that process, the transverse shear stress x, should be 

evaluated from a trial value of X by assuming that the Ax arises from wholly- 

elastic behavior of Ay, but Ay is always zero and hence X is also zero). To 

avoid an unduly complicated analysis, and still achieve good accuracy, and also 

because the transverse shear stress, X , often may be small compared with the 

axial stress a, (Ref. 20) , the following yield function is adopted for the 

present Bernoulli-Euler-type beam 

2 

- a ~ 

However, the yield function (Eq. 4.42) is used for the Timoshenko- type beam 
where both the axial and transverse shear stresses and strains are taken into 
account; this will be described in the next subsection. 



An illustration of the method of computing the axial stress and/or 

plastic strain increment is presented as follows. One begins by knowing the 

sublayer stress 0^ at time t^_^ for the kth sublayer of the jth depthwise 

Gaussian station, and the strain increment Ae. . at station j at time t. 

3 i 


(that is, the strain increment from time t i _ 1 to time tj 


One then takes a 
which is computed by assuming an elastic 


trial value (superscript T ) of a.. . 

+ 3k, i 

path ; 

= V*'-' +E "V 

A check is then performed to see what the correct value of 0 . 

3k, i 


(4.44) 
must be. 




then 


= °jh 

and 

o 

II 

. i 

*< 

<1 

If 

^ ^ok 

then 


= r 0 k 

and 

7 M £ 

If 


then 

°fk.i 


and 

E 


(4.45) 


where E is Young's modulus. 


This procedure is applied to all sublayers of each Gaussian station j ; 

+ It should be noted that the subscripts in quantities such as ^ ^ , for example , 
represent only identifying labels, not tensor notation as used heretofore. 
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having done this, the axial force and moment of the beam cross section can be 
determined by 


L = / ( <r d A = b J XL ( 




A f k > 


M = f(a^dA-bf CjkAn) 

A 1 * « 


(4.46) 


In a similar manner the integration of the plastic strain over the cross 
section of the beam element can be determined by 


!\ 4 2 (E. t' k Aft ) 

5,-(£ £/„ A,, ) 


(4.47) 


where b is the width and h is the thickness of the beam and A is a corabina- 
tion of the mechanical sublayer weighting factor and the Gaussian weighting 
factor W . , which is defined by 

( E k - E 


fc- 


(4.48) 


In Eq. 4.48, W. is the Gaussian weighting factor (Ref. 131) and 

■ 3 r- = - £k 

t ’ k t'k - £ /f-> 

is the kth slope of the polygonal approximate stress-strain diagram, 
be verified that the relations . 


(4.49) 
It can 


( 


are satisfied. 


4-V/H = 


"* vv j-' r i' An o~ k 


If desired, the sublayer yield stresses may be treated as strain-rate 
dependent. Since the strain increment at the jth Gaussian station and hence 
the strain rate is known at this stage of computation, then the rate-dependent 
yield stress a ^ of this kth sublayer at station j is 


fy “ ^ok [ 


b£/it | f 


l 

r ] 


(4.50) 


where D and p are empirically-determined constants for the material 
and may, in general, be different for each sublayer. 


ok 


is the static uniaxial yield stress of the kth sublayer 
at any jth Gaussian station 
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4.4 The Timoshenko-Type Curved Beam Element 


4.4.1 Displacement Field 

In the previous subsection, the derivation of the beam element properties 
is based on the Bernoulli-Euler hypothesis in which the transverse shear deforma- 
tion is assumed to be zero. In the present subsection, the formulation for a 
general curved beam element with nonzero transverse shear deformation taken 
into account (termed a Timoshenko-type beam element) will be presented. How- 
ever, the following assumptions (Ref. 1S9) are used: plane cross sections per- 
pendicular to the undeformed centroidal axis remain plane and suffer no strain 
in their plane, although they no longer remain perpendicular to the deformed 
centroidal axis. 


With these assumptions, the displacement field, including transverse 
shear deformation of the beam may be specified by the middle plane displace- 
ments and cross-section rotation, as follows: 


7 i?, s ) = vi*i) + 5 e 
w (7 , 5 ; = wc^ > 


where v(n) 
w(n) 
6(n) 
C 


3 axial displacement of the middle plane 
= transverse displacement 

= rotation of the plane cross section about the x axis 
= normal distance from the centroidal axis (middle plane) 


4.4.2 Strain-Displacement Relations 


Neglecting the variation of the transverse shear strain across the 
thickness of the beam, the expression for the engineering components of the 
strain distribution may be written as 


£ ( S ) = £. + C } 

><>?, C ) = ~Yo < ? > 


(4.52) 


where is the middle-plane axial strain 

k is the curvature change 

Y is the transverse shear strain 
'o 
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The strain-displacement relations of the curved-beam element may be 


expressed as 

e. 


, XL 
?? R 



K 


Y. 


9 9 

*7 


= ( 


l tv _ X. 

R > 


■+ e 


(4.53) 


In this equation also, the only nonlinear term retained is due to the rotation 
of the centroidal axis about the x axis. Equation 4.53 may be expressed in 
matrix form as 


2. 

3 

- L-37 
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i v l 

0j in 
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1 ( R 
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*7 
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1 - f 
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O 
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d 

37 

1 _J j 

1 w 
1 6 

—1 

CQ 

_J 

III 

1 

J 


(4.53a) 


It should be noted that in the previous Bernoulli-Euler-type beam ele- 
ment, the highest derivatives upon which the strain depends are the second 
derivative of the transverse displacement, w, and the first derivative of in- 
plane axial displacement, v. In order that the assumed displacement inter- 
polation function be admissible in the variational argument of Eq. 3.3, it is 
required that the assumed displacement function of w at least possess a second 
derivative and v possess a first derivative; this means that the assumed dis- 
placement field must generate continuous displacements and continuous normal 
slopes at the interelement nodes. However, the inclusion of transverse shear 
deformation reduces the order of the derivative requirement and hence also the 
the stringency of the compatibility imposed on the assumed displacement field, 
because the strains depend only on the first derivative of displacements v,w 
and the rotation 6. For the variational argument of Eq. 3.3 to be defined. 
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it is required that the assumed displacement field have at least a first deriva- 
tive and provide continuous displacements and rotations at the interelement 
nodes. Consequently, the simplest assumed displacement field for this Timoshenko- 
type beam element with the rigid-body displacement modes included is 


v' 


cos 4> 

SIN <j> 

- zcos (p-tysiN <£ 

y o o v 


A 1 

o 

< w 

- 

~SIN $ 

cos4> 

z sin4 +y cos 4> 

0 0 0 

< 

Pi 

■ 0] 


0 

0 

0 

o / 7 , 


,k 



(4.54) 

The generalized nodal displacements {q} are defined to b <i the three degrees of 
freedom v, w, and 6 at each node of the element as follows: 

\ - L V V ^ V A+I J (4.55) 

In a manner similar to that described in the previous subsection, one 
may write 

f ft j = C A J {p} (4.56) 


and 



= [N<7>] ffj 


(4.57) 


It perhaps should be mentioned tliat by this linear interpolation function 
(Eq. 4.54) of the displacement, the strain and moment representation will be 
very crude unless the element size is kept small enough, since as can be seen 
in this formulation, the bending strain is constant over each individual ele- 
ment. In order to improve the strain representation, higher order displace- 
ment interpolation functions and hence more degree of freedom (or internal 
nodes) should be used (see Ref. 164) . Further discussion of this matter is 
given in Subsection 5.2.2. However, for the purpose of illustrating the cal- 
culation procedure, the linear interpolation function is discussed here. 
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Combining Eg. 4.57 with Eq. 4.53a results in 

e . = ld,j tp * j - Lpm\iv tl \p 

K = lD 3 J \f) 

~Yc = if} 

where 

L^/J = lB/j E N C ? J J i=\.Z.3,4 


(4.58) 


(4.58a) 


Also, the strain increments from time t^ 


to time t. are given by 


where 

— T 

^■= -I ?*£}*• 
= lAj 


4,4.3 The Discrete-Element Property Matrices 

The element mass matrix and the element generalized nodal forces can 
be obtained following the same argument as in the previous subsection. Also, 
in addition to L and M, the transverse shear force, S, of the cross section, 
A, is given by 

c = ((-i d A 

° J ) (4.60 

A 

where T is the transverse shear stress. 


Then, from the Principle of Virtual Displacements equation, Eq. 3.14, 
one may obtain for the improved formulation : 

Z ( OJ \ 'j\ + ff S + [M {fj - if 1 ) = o ( 4.6 
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where 


I > 1 

If) 
I PI 

[hi 



= /^'(|D,)L * ( Dj i M -^fD4S )<<? 


(4.61a) 


Also, for the conventional formulation : 

t L5p([mlfj] - rklffrt - if)- 

where 

ft J =f !i "( jx>,j eU lAj + fa j ^ J L qj - [qj & U lQj ) ^ 
H>} = // /^'(Ee^hl-CE^J+s^fn}) dU„ 

/l A 

4 A' 


If £) = £fcf> <hp f D 4 i t \ ) (d,} ^ 

1^ 


Finally, one can recast Egs. 4.61 and 4.62 in terms of global generalized 
displacements {q* } to obtain the corresponding equations of motion as described 
in Subsection 3.2. 
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4.4.4 S tress-Strain Relations 


For this Timoshenko-type beam element with transverse shear deformation 
included, the calculation of stresses and stress resultants for both normal 
and transverse shear components is illustrated as follows. 

One begins by knowing the sublayer normal stress 0^ and transverse 
shear stress at time for the kth sublayer of the jth depthwise 

Gaussian station (layer) and the jth layer normal and transverse shear strain 
increments A E. . and Ay. . , respectively, at time t. . The trial stresses are 
calculated from the relations 




= * E 

-V'- ^ 


(4.63) 


Then, the trial stresses are introduced into the Mises-Hencky yield function 




(4.64) 


where E is Young's modulus of elasticity 
G is shear modulus of elasticity 
and 0 Qk is the yield stress of the idealized elastic perfectly- 
plastic kth sublayer. 

T 

If 4^ ^ 0, the stress state lies within the yield surface, no plastic 
flow occurs within this time step, and the actual stress increments arise from 


wholly elastic behavior, 

then 




®l'k, * 

a jk,. 

"T., • 

’ 1 k,A 

= Tr.T • 

(4.65) 

and 





■ = 


> -rju 

= -y> 

(4.66) 

T 

If, on the other hand 4\ 

> 0, plastic yielding has 

occurred. Then, 

the flow 


rule. Eg. 2.82a gives the plastic strain increments as 
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a e: 


yk , / ^ ^ ^y/c, / 




(4.67) 


and the actual stresses are 


&jk,i - ( !~^E) cr- fc v , ^,, =(,_3A (4-63) 

where X* = 2A and can be solved from the requirement that the actual stress 
state must be on the yield surface. Thus, the following condition must be 
satisfied: 


* =f< * 3 - a;* 2 ) =0 

Substituting Eq. 4.68 into Eq. 4.69, one may solve for X* as follows: 


A* = 


where 


B + 

a = [e ’(<£.,)*- 

B - [E ((£.,)' * 1 


(4.69) 


(4.70) 


C = 


(4.70a) 


With X* obtained, the stress state at time t^ and the plastic strain increment 
from time to t^ are known. This process must be carried out for each 

layer (i.e., depthwise Gaussian station and sublayer). Once the stresses in 
each layer and sublayer have been determined, the axial force, moment , and trans- 
verse shear force of the cross section can be obtained by 
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and the integration of the plastic strain over the cross section can be ob- 
tained by 

(j £*dA ef k A fk ) 

(( st**A S, 

(( ytjA - kf r (Z-y* k Af k ) ,4., 2 : 

A 1 K 

Then the equilibrium equation, Eq. 4.61 or Eq. 4.62, must be used next in the 
calculation cycle to find the displacement or displacement increment of the 
next time cycle, as discussed in Subsection 3.3.1. 
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SECTION 5 


EVALUATION AND DISCUSSION 


5.1 Introduction 

In Section 4, the curved beam element, which may undergo large deflec- 
tions and elastic-plastic strains, as well as deforming such as either to in- 
clude or omit transverse shear deformation, has been developed for the dis- 
placement variational model. The timewise numerical 3-point central-difference 
finite-difference procedure is employed to solve the resulting system of coupled 
second-order ordinary differential equations. 

In order to evaluate the accuracy and versatility of the present finite- 
element formulation and solution scheme, this analysis has been implemented in 
a computer program and several numerical examples have been carried out. First, 
in Subsection 5.2, comparisons are made between the present finite-element solu- 
tions and known analytical solutions for small-deflection linear-elastic tran- 
sient responses of mechanically- loaded beams. Next, in Subsection 5.3, the 
present predictions are compared with those from available finite-difference 
(both spatial and temporal) predictions and with experimental observations for 
large-deflection, elastic-plastic transient responses of impulsively- loaded 
beam and ring structures , and various features of the present method are 
assessed. 

5.2 Small-Deflection Linear-Elastic Transient Responses 
of Mechanically Loaded Beams 

In order to check on the proper functioning and correctness of the 
present analysis and computer program, the small-deflection linear-elastic 
transient responses of beams have been analyzed first; the finite-element pre- 
dictions have been compared with available analytical solutions. Two beam 
problems have been studied: one pertains to Bernoulli-Euler (or Kirchhoff) 
deformation behavior while the other includes a significant amount of transverse 
shear deformation. These examples are discussed in the following. 
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5.2.1 Small-Deflection Linear-Elastic Responses of a Beam 
with Bernoulli -Euler-Type Deformation Behavior 

A simply-supported beam is subjected to a uniform lateral transient 
loading over the whole spam . The geometry , dimensions , material properties , 
and transient-loading history are depicted in Fig. 6. The Bernoulli-Euler- 
type small-deflection elastic beam element with a linear interpolation function 
for the axial displacement v and a cubic interpolation function for the trans- 
verse displacement, w, was employed; this element has been termed an LC element 
in Section 4. Because of symmetry, only one-half of the span of the beam was 
modeled; five equal- length discrete elements were used in the attendamt finite- 
element analysis. The solution was obtained by using the 3-point central -differ- 
ence timewise integration method with a time increment size of At « 4 usee which 
satisfies the stability criterion. At <_ (2/“^) , for this method (see Ref. 148) 
where tii represents the largest natural frequency contained in the mathematical 
model, (M]{q*} + [K]{q*} = 0 which approximates the actual linear elastic small 
deflection structure. 

A comparison of the mid-span transverse deflection response predicted by 
using the present finite-element scheme with the exact normal-mode solution is 
shown in Fig. 6. It is seen that very good agreement between these solutions 
is observed. It should be noted that for this small-deflection linear-elastic 
straignt beam with Bernoulli-Euler deformation behavior: (1) the inplane 
(axial) displacement is zero and (2) the selected assumed cubic displacement 
function for the transverse displacement w is, in fact, identical with the / 

exact displacement field. Hence, the finite-element calculation which utilizes 
the central-difference time integration method gives very accurate amplitude 
and phase predictions for small-deflection linear-elastic behavior as long as 
the time increment size used is small enough to satisfy the stability criterion. 

5.2.2 Small-Deflection Linear-Elastic Transient Responses 
of a Beam with Timoshenko-Type Deformation Behavior 

The second example is selected to test the convergence of predictions 
utilizing the various assumed displacement functions for the Timoshenko-type 
beam element when it is applied to a small-deflection linear-elastic dynamic 
system. Transverse shear deformation and rotatory inertia effects are included 
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in the formulation. A free-free beam is subjected to an applied loading con- 
centrated at the mid-span with a triangular pulse time history (this problem 
has been analyzed in Ref. 165) . Figure 7 gives the geometry and the dimension- 
less parameters for this problem. 

Because of symmetry, only a half-span of the beam was treated in the 
finite-element solution. The Newmark constant average acceleration timewise 
integration scheme + (Y = 1/2, 3 = 1/4) with a dimensionless time increment 
size At * 0.0025 was used, where At - At J El/pbh /£ . It should be recalled 
that this integration operator is unconditionally stable for linear (small- 
deflection) elastic dynamic systems; however, too large a At may cause some 
phase shift from the correct behavior. 


For this beam problem with a significant amount of transverse shear de- 
formation, the following four types of assumed displacement fields (designated 
as T1 through T4) have been tested (note that zero inplane displacements v are 
involved — v is ignored) : 


(Tl) Linear functions in S' for both transverse displacement w, and 
rotation 6: 

W = 0., -+ !Xi S, 

e - a 3 e, 

The generalized coordinates {q} are selected such that there 
are two degrees of freedom (w, 0) at each end node i and i+1 
of the element: 

\f\ - L QjC • V/*V, 0/*. J T 

(T2) Cubic variation of w and linear variation of transverse 
shear strain 

w = a, a t H, f a s x* -+ a 4 ^ 


(5.1a) 


(5.1b) 


(5.2a) 


This operator rather than the 3-point central-difference operator was employed 
in an attempt to use a larger At than the latter permits, and thus reduce the 
computing time. 
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(T3) 


and 


(T4) 


Three degrees of freedom (w, , 8) are selected at each end 

node of the element and the generalized coordinates are : 








0 / V I _J 


The reason for choosing the linear function for the transverse 
shearing strain (i.e., quadratic function of the rotations) is 
that the bending effect dominates the transverse shearing ef- 
fect when the element size is large, and the bending strains 
are derivatives of the rotation. Accordingly, this function 
could represent the strain energy of the element more accurately 
for a large mesh size (Kef. 164) . 

The same assumed displacement functions as for T2 (i.e., with a 
cubic variation of w and a linear variation of except that 

the generalized coordinates are selected such that there are 
two degrees of freedom (w, 8) at the two end nodes and at a 
midpoint node of the element. Thus, 

.2 . „ *3 


w = a, + a i £ 


a. 


a* 




% 


(5.2b) 


(5.3a) 


W. 


\f\ = l w* Qx ^ e m 

Quadratic variation of w and linear variation of y 


w = a, 
^ 5 " ° A 


a, i + a 3 V 




w 


Two degrees of freedom (w, 8) at the two end nodes and one 
degree of freedom, w, at the midpoint node of the element are 
selected. The generalized coordinates are 


lfj = L K- 8 .- 




XM 


irl 


e A^I J 


(5.3b) 


(5.4a) 


(5.4b) 
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Shown in Fig. 8 are the corresponding dimensionless transverse shear 
force responses preuicted at tne quarter-span station by using the four differ- 
ent sets of assumed displacement functions iTl tnrougn T4) with various numbers 
of discrete elements for the naif span. Shown also are the "exact" modal solu- 
tions of Ref. 165 based on the Timoshenko assumption for convenient comparison. 
It is seen that the use of linear interpolation functions for both w and d (Tl) 
gives very crude coarse-mesh transverse shear force responses compared with that 
from the modal solution. Whereas, predictions obtained by using the nigher- 
order interpolation function (T2, T3, and T4) with a coarse-mesh finite- 
element array show better agreement with the modal solution than the predictions 
by using linear interpolation function (Tl). However, as the number 
of finite elements employed increases (the element size decreases) , all four of 
these interpolation-function predictions converge to the modal solution both 
in pnase and in amplitude. If one bases the comparison on the total number of 
degrees of freedom (unknowns) which were used in the finite-element solutions, 
it is seen that, for a given number of unknowns, the predictions obtained by 
using T2 or T4 type interpolation-function elements leads to a solution which 
is closer to the modal solution than is the case if Tl or T3 type interpolation 
function elements are employed. The reason for this comparative behavior is 
that in the Tl type (linear) interpolation function, the strain and moment 
representation over each element are very crude; whereas in T2, T3, and T4 type 
interpolation function elements , the strain and moment representation is much 
improved over that with the Tl type element. Also, by using T3 type elements, 
the mesh size is relatively larger compared with that from using T2 or T4 type 
elements if they have tne same number of degrees of freedom for the half-span 
of the beam. Finally, it should be noted that the size of the finite elements 
which provide good shear response agreement with the exact solution are such 

that their (equal) length is less than the depth of the beam (see Fig. 3) 

such a required modeling pertains to problems which include transverse shear 
and rotary inertia. 

It perhaps should be noted that under various loading conditions, slope 
discontinuities along interelement boundaries are permitted when considering 
the presence of transverse shear deformations ; thus , the T2 type element 
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over-constrains tne continuity of slope at the interelement nodes. Based on 
this consideration, it may be concluded that the T4 type assumed-displacement 
function is the most efficient one and the T2 type assumed displacement func- 
tion coraes next, when they are applied to beams with transverse shear deforma- 
tion benavior. 

The finite-element predictions for the quarter-span moment responses and 
tne midspan deflection responses are shown in Figs. 9 and 10, respectively. 

The above-mentioned convergence behavior of the shear-force responses obtained 
by employing the four types of assumed displacement functions are also observed 
in both of these latter two figures for the moment and displacement responses. 

This example illustrates that the linearly-varying (Tl) assumed displace- 
ment Timoshenko- type beam finite element can provide accurate transient response 
predictions only if the element size is kept small enough. However, in order 
to obtain more accurate coarse-mesh solutions, one would need (a) to employ 
iiigner-order assumed displacement functions (T2, T3, T4) or (b) to use an 
assumed stress hybrid finite-element model (Ref. 166) . 

5.3 Large-Deflection bias tic-Plas tic Transient Responses 
of Impulsively-Loaded Simple Structures 

In order to evaluate the reliability and accuracy of the present finite- 
element metnod for predicting large-deflection elastic-plastic transient re- 
sponses of simple structures, the various features and options of the present 
prediction metnod are examined in this subsection. Also, comparisons of the 
present predictions with other available finite-difference (both spatial and 
temporal) predictions (Ref. 44) and with experimental observations (Refs. 44, 

167 and 16U) are made. 

5.3.1 Example Problems Analyzed 

As examples with which the various features of the present finite- 
element prediction method can best be illustrated, the following three types of 
problems have been analyzed; the geometries and the types of finite elements 
employed are presented in Fig. 11 for reference convenience. 
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(1) A straight beam of 60G1-T6 aluminum alloy is clamped at 
each end and has dimensions: 1/8-in. thickness, 1.2-in. 
width and 10-in. span between supports . The beam is 
loaded impulsively over a spanwise segment centered at 
the midspan and covering a length of 2 inches, as de- 
picted in Fig. 11a. 

(2) A 6061-T6 aluminum alloy ring sector with a subtended 
angle of 315° is clamped at each end and has dimensions: 

2.935 in. mean radius, 0.123-in. thickness, 1.197-in. 
width. The clamped ring is loaded explosively over a 
peripheral sector of 120° of its exterior, as shown in 
Fig. lib. 

(3) A free circular 6061-T6 aluminum alloy ring, 0.124-in. 
thick, 1.195-in. wide, and 2.937 in. mean radius is 
subjected to severe explosive loading over a peripheral 
sector of 120° of its exterior, as shown in Fig. 11c. 

As a matter of convenience for reducing the computer time and storage 
required by the solution of these problems, these three examples were treated 
as symmetrical problems. Taking account of the symmetry of the impulsive 
loading, geometry, and boundary conditions, only half of each configuration 
was modeled in all of the attendent discrete-element analyses. 

5.3.2 Effects of Using Various Timewise Integration 

Operators: Central-Difference Method, Houbolt's 
Method, and Newmark's 3-Method 

In this subsection, the numerical stability behavior of the timewise 
integration operators: central-difference method, Houbolt's method and Newmark's 
p-method (y = 1/2, p = 1/4) when employed for the solution of large-deflection 
elastic-plastic transient responses of tne impulsively loaded clamped beam will 
be studied. The conventional finite-element formulation of the equilibrium 
equations is used, where the large-deflection and elastic-plastic effects are 
accounted for through the use of "equivalent generalized forces" which are given 
automatically from the variational formulations. 
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To minimize the computer time for the structural response calculations, 
one should use the largest permissible time increment which will avoid numeri- 
cal instability (for example, roundoff error blowup or truncation error accumu- 
lation) and still provide a reliable prediction. Unfortunately, for the 
present nonlinear dynamic system, a reliable and validated criterion by which 
tiie proper time-step size can be chosen a priori is not readily available for 
any of these methods. 


If the 3-point central-difference (timewise) method is used, as was 
pointed out in Subsection 3.3.4, the judicious selection of the proper time 
increment At can be guided by the stability criterion of a corresponding linear 
dynamic system, At <_ , as an initial selection ; numerical experimenta- 

tion then subsequently can provide the suitably smaller At to insure stability 
where ^ max represents the largest natural frequency in the mathematical model, 
[M]{q*} + [K] {q*} = 0, which approximates the actual (linear-elastic small- 

deflection) structure. Thus, it would be very valuable to know , such that 

max 

the "suitable initial At" can be chosen immediately and hence reduce the amount 
of subsequent numerical experimentation. Figure 12a presents of the 

clamped beam as a function of the number of elements per half span; the 
Bernoulli-Euler-type LC beam element with a consistent mass matrix is used to 
model the structure. The maximum u was obtained by an iteration process in 
double precision applied to (see Art. 4.5 of Ref. 27) : 


))f! (5.5 

Also shown in Fig. 12a are the maximum frequencies of pure membrane behavior 
and pure bending behavior (with or without including the rotatory inertia 
effect) . It is seen that the maximum frequencies of combined membrane and 
bending behavior are equal to the pure membrane maximum frequencies, when the 
mesh size £, to thickness, h, ratio is large (£/h _> 4), and are equal to the 
pure bending maximum frequencies when the mesh size to thickness ratio is 
small (e/h 4) . This is to be expected for the linear Bernoulli -Kirchhoff 
beam system, because membrane and bending effects are decoupled. It is also 
noticed that the rotatory inertia effect arising from the consistent mass 
matrix can be significant when the mesh size is small, and the neglect of 
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rotatory inertia by deleting those terms from the consistent mass matrix leads 
to a maximum frequency higher than that obtained by including this rotatory 
inertia effect (as seen in Fig. 12a) . 


It should be noted that simple analytical methods to estimate the upper 


bound and the lower bound of ^ aax have been presented, for example in Ref. 162, 
in terms of the maximum eigenvalue of all the individual element matrices and 
its associated eigenvector: 

s a/ . * t > 

(5.6) 


< U) „ 
__ / •* r \ max 

+ r l* j [m n : |xj l 


2 2 2 
where max(v ) denotes the largest \> of all the finite element and V is the 
n * n n 

maximum eigenvalue of the equation 


p z fmj = [IcJ !*} 


(5.6a) 


and 


{x*} is the eigenvector of Eq. 5.6a for the element j, 
2 

which has max (v ) , 
n 

£' indicates summation over only the neighboring 
elements of element j. 


The element matrices [m ] and [k ] involve fewer degrees of freedom than 

n n 

the assembled matrix [M] and [K] , so it is relatively much easier to find the 

element's maximum eigenvalue than the co of the assembled matrix representing 

max 

the complete structure. However, the bounds may not be very sharp; also the 

boundary conditions (which the cited "bound method" does not take into account) 

and the nature of the problem will affect the w as shown in Fig. 12b where 

max 

the upper and lower bounds of the maximum frequencies of pure beam bending 
behavior (by using consistent mass matrices from which are deleted the rotatory 
inertia effect) are given by + (Ref. 162) : 


900 


EI 


$ 


a), 


max 


^ 8 A-00 p g4 


(5.6b) 


where p is the mass per unit length and £ is the element length. In view of 


Similar bounds could be developed when one uses lumped mass matrices such as 
those discussed in Subsection 4.3.5. 
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the wide bounds given by these results and in order to reduce the effort in- 
volved in subsequent numerical experimentation to determine a maximum permissible 
At, the exact w max has been evaluated in the present analysis; fortunately, 
for the beam (or ring) problem, the total degrees of freedom are not too many 
to be handled by the scheme represented by Eq. 5.5. However, for problems with 
an enormous number of degrees of freedom, one would perhaps need to resort to 
the simple analytical bound methods to estimate the upper and lower bounds of 

0 ) ; for such cases, the “bound scheme” may be more economical and efficient, 

max 

Now, turning to the large-deflection elastic-plastic transient response 
predictions, the half-span of the beam is modeled by 10 elements and the time- 
wise 3-point central-difference operator is used. The critical At, if based 
on the stability criterion of a corresponding linear system would be 
At „ = 2 /oj = 1.47 Msec. Computational experiments have been carried out 

CX, IT13.X 

using various time step sizes as shown in Fig. 13a to predict the midspan de- 
flection responses. This clearly demonstrates the immediate divergence of the 
predictions if At is only slightly greater than At^tAt = 1.5 Msec = 1.02 x At^) ; 
reliable predictions are obtained if At <_ 0.99 x At^ = 1.45 Msec. Calculations 
also have been carried out by using Houbolt's method and Newmark's method as 
shown in Fig. 13b and 13c, respectively; it is observed for both of these 
methods that for At values which are too large , the predicted response degrades 
gradually but badly from the correct behavior. The critical At for reliable 
predictions was found to lie between 6 and 8 Msec for Houbolt’s method and be- 
tween 3 and 4 Msec for Newmark's method. 

In view of these results and those of Refs. 55, 75, and 156, it is con- 
ceivable that the following situation may generally be considered to be true. 

The introduction of material nonlinearity often decreases the highest natural 
frequency of the system because the plastic pulse travels at a velocity which 
is less than the elastic pulse velocity, but elastic response contributions 
are still present and govern the allowable At. However, the geometric non- 
linearity effect (large deflections) renders Houbolt's method and Newmark’s 
method no longer to be "unconditionally stable". 
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It should be remembered that for this beam problem the calculation of 
strain increments from displacement increments and displacements, Eq. 4.25, at 
any time t is given by 


(A£) = ( A £.) + ^ 


where 


(A £ ) 

v O 'tn 


- tl±l) +12*') (2 *2) --Lfi**') 

' 3 73 'm K d <2 'm ' d '> *1 2 V 3- 7? A 




L (2* vs * 1 
2 ( Fy - m 


K = - ( 


_ / ^_a_vy J 

'» 


(5.7) 


2 

The higher-order term l/2(3Aw/3^) m has been included in the above calculations. 
The predictions made by neglecting this term in the calculation of the strain 
increment are shown in Figs. 14a, 14b, and 14c for the central-difference me- 
thod, Houbolt's method, and Hewmark's method, respectively. Comparing Fig, 14 
with Fig. 13, it is seen that the neglecting of this higher-order term may de- 
grade the long-time responses due to the accumulation of "errors of approxi- 
mation" introduced at each time step, especially for larger At as can be seen 
more prominently in the predictions obtained by Houbolt's method and Hewmark's 
method. Accordingly, it is recommended that the exact strain-increment equa- 
tion (Eq. 4.25) including all the linear and nonlinear terms in the displace- 
ments and the displacement increments should be used; fortunately 
the computer time and storage increase is insignificant. 


It should be noted that Houbolt's method and Hewmark's method are im- 
plicit in nature; that is, the generalized nodal forces (which may be due to 
large-deflections and elastic-plastic effects) at each time step depend on the 
displacements (or stress, strain) at that time step, which remain to be de- 
termined; thus, iteration or extrapolation is needed at each time step. Linear 
extrapolation (see Ref. 75, for example) by using the generalized nodal forces 
at two previous time steps is employed in the present calculations. The central- 
difference method on the other hand is explicit in nature and thus no iteration 
or extrapolation is required at each time step. The storage of the displacements 
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at three previous time steps is required by Houbolt's method but information 
only at two previous time steps is needed when the 3-point central-difference 
method or Newmark's method is used. Also, for At values which are too large, 
the very gradual degradation of the responses predicted by using Houbolt's 
method or Newmark's method gives no warning to the analyst that this degrada- 
tion may be happening. However, the central-difference method usually exhibits 
a violent degradation of response when At is too large, thus warning the analyst 
of this state of affairs. 

Based on these considerations*, only the timewise central-difference 
method is employed in the following example calculations. However, it should 
be mentioned that based on the present information, it is still far from con- 
clusive as to whether any one timewise operator is superior to the others for 
analyzing nonlinear transient response problems of the present type. 

5.3.3 Comparison of the Conventional Formulation Versus the 

Improved Formulation for the Dynamic Equilibrium Equations 

By using the timewise 3-point central-difference method, comparisons 
have been made of the responses of the impulsively- loaded clamped beam ob- 
tained by employing the conventional finite-element formulation with those ob- 
tained by using the improved finite element formulation of dynamic equilibrium. 
Complete agreement of these predictions is observed in Fig. 15 for both the 
large-deflection elastic-plastic transient responses and the small deflection 
linear-elastic transient responses. However, as was discussed in Subsection 3.3, 
the improved formulation shows significant simplification in form over the con- 
ventional formulation for solving large-deflection elastic-plastic dynamic equi- 
librium behavior. Also, the computer storage and manipulations required for the 
improved formulation are less than those required for the conventional formula- 
tion. For this problem, it has been found that when using the conventional 
formulation, the computer time is about 24% more than that required when the 


+ Also in order to make convenient comparisons, because in the available inde- 
pendent finite-difference (both spatial and temporal) predictions involving 
large-deflection elastic-plastic behavior (Ref. 44), the timewise central- 
difference method is used. 
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improved formulation is used. Based on these considerations, only the improved 
formulation will be used in the following further example calculations of large- 
deflection elastic-plastic structural transient responses. 

5.3.4 Comparison of LC Versus CC Assumed Displacement 

Functions for the Bernoulli -Euler-Type of Ring Element 

The free ring has been analyzed by using the Bernoulli-Euler type of 
ring element with either CC or LC assumed-displacement functions. Shown in 
Fig. 16a is the maximum natural frequency to for the linear behavior of the 
finite-element representation of the ring as a function of the number of ele- 
ments per half ring. It is seen that for the same mesh size, the use of CC 

assumed-displacement elements has a larger u) (hence requires a small critical 

max 

time step At 0 = 2/ui > than that occurring when the LC assumed displacement 
cx max 

elements are used. The ring centerline separation time histories predicted by 
using CC assumed-displacement elements compared with that predicted by using 
LC assumed displacement elements are shown in Fig. 16b, where the ring material 
is considered to be elastic linear-strain hardening (EL-SH) . The experimentally 
observed response is also shown for convenient comparison. It is seen that as 
the structure is modeled as more and more finely subdivided, the LC element so- 
lutions converge and provide a somewhat stiffer response compared with experi- 
ment, while the solutions obtained from using the CC element converge more 
rapidly but tend to be "too flexible". However, the strain- time histories as 
shown in Fig. 16c indicate that the strain responses predicted by the CC ele- 
ment are very close to measured values , while the LC element under-predicts 
the strain. 

It should be noted that for this free ring subjected to severe impulse 
loading (see Fig. 11c) , one may expect the strain-rate effect to become rather 
important. The central line separation and strain time histories predicted 
from the elastic linear strain-hardening and strain rate (EL-SH-SR) calculations 
are shown in Fig. 17a and Fig. 17b, respectively. Far better (in fact, excellent) 
agreement with experiment of the CC element predictions than those of the LC ele- 
ment predictions are observed. It should be mentioned that for this free ring 
involving a rather severe degree of response, computational experiments have 
indicated that reliable deflection and strain predictions are obtained if the 
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time steps are At < 0.8 x At „ = 0. 8 ( 2 / 0 ) ). Comparing with At < 0.99{2/u) ) 

— ex. max — max 

for the clamped beam, it is believed that the critical time step is affected 
by the severity of the structural responses. That is, large deflections play 
the key role in effectively stiffening the structure and thus requiring a 
smaller At to avoid round-off error instability. 

The above-mentioned comparable behavior of the CC element predictions 
versus LC element predictions are also observed in the clamped ring calculations; 
see for example, the central deflection responses presented in Fig. 18. 

In view of the present results and those of (Ref. 161) , it can be con- 
cluded that the CC-type assumed displacement function exhibits significantly 
improved predictions , especially for the strains , over the LC assumed-displace- 
ment elements. Also, the former converges very rapidly but at the expense of 
a smaller allowable time-increment step compared with the latter. 

It should be noted that in the above example, the converged solution 
obtained by using a finer mesh size is caused not only by the better . approxima- 
tion of the original continuum, but also by the better representation of the 
abrupt change of the initial impulse loading imparted to the system at the edge 
of the high explosive. This matter is discussed further in Subsection 5. 3. 6.1. 

5.3.5 Comparison of the Use of a Consistent Mass Matrix 
Versus a Lumped Mass Matrix 

If the lumped mass matrix is used (a = 1/2, 5 = 1/24, see Eq. 4.32) for 
the analysis of the clamped beam example, it is observed in Fig. 19a that, for 
the same mesh size , the maximum frequency represented by the lumped mass matrix 
system (linear-elastic, small deflection) is smaller than that obtained by the 
use of a consistent mass matrix system. Hence, a larger time increment size 
(to avoid numerical instability) can be used for the lumped-mass-mat'rix system 
than for the consistent-mass -matrix system. 

Figure 19b shows the midspan deflection responses for large-deflection 

elastic-plastic strain-rate dependent behavior. The responses predicted by 

both types of mass-matrix systems are quite close to each other, where the 

half-span of the beam is modeled by 10 elements and the time increment size 

used for a stable solution is At < 0.99 At . = 0.99(2/u ). That is, 1.45 Msec 

— c x. max 
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for the consistent mass matrix system and 2.45 Msec for the lumped mass matrix 
system. A violent degradation of the responses occurs when At is only a trifle 
larger than At Comparisons of lumped mass model (a = 1/2, 6 = 1/24, see 
Eq. 4.32) predictions versus consistent mass model predictions for the free 
ring example using Bernoulli-Euler CC-type ring elements are presented in 
Figs. 20a, 20b, and 20c for the maximum natural frequencies, the central-line 
separation, and the strain responses, respectively. Again, by using the 
lumped mass matrix, a smaller maximum natural frequency and good accuracy of 
the predicted nonlinear responses compared with the use of the consistent mass 
matrix are observed. 

It should be noted that the use of the lumped mass matrix (with a = 1/2, 

6 = 1/24) not only decreases the maximum natural frequency (hence enables one 
to employ a larger time increment step for the response calculation than the 
consistent mass matrix system) , but also reduces the storage and computer time 
required for the solution of the transient response problem. Because the lumped 

mass matrix is a diagonal matrix, [M] = [m.6. .], its inverse is just 

—1 ^ 

(M] => [1/nu ) ] . However, it should be noted that further studies need to be 

conducted to develop mass matrix lumping rules which are appropriate for various 

user-selected criteria. 

5.3.6 Assessment of Some Features of the Method 

Among the various features of the present finite-element analysis which 
are examined and discussed in the following are: 

(a) the effects of various initial velocities, specified 
at the nodal points of the finite-element assembled 
structure, to approximate the impulse loading which 
is produced by the detonation of a sheet of finite- 
span high explosive; 

(b) effects of the number of spanwise Gaussian points 
used to evaluate the properties of each discrete 
element {p} and [h] , and the number of aepthwise 
Gaussian points used to evaluate stress resultants 
(axial force, moment and/or shear force) at each span- 
wise Gaussian station; 
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(c) tiie effect on the predicted transient response of 
including strain-rate sensitivity; 

(a) the effect on the predicted transient response of 
including transverse shear deformation. 

Each of these matters is discussed as* follows . 

5. 3. 6.1 Effects of Using Various Initial Nodal Velocities 
to Approximate the Impulse Loading 

Experimentally, the impulsive loading may be produced by the detonation 
of a sheet of high explosive (HE) . Between the HE and the test specimen is a 
thin layer of a suitable buffer material to prevent intense stress-wave- 
induced spall fracture of the test specimen. Experiment indicates that a 
nearly uniform initial normal impulse is imparted to those portions of the 
specimen immediately underneath the HE layer. However, for the region of the 
beam near the spanwise edges of the HE layer, a very steep gradient of imparted 
impulse is observed. A typical normalized distribution (Ref. 44) of the im- 
parted impulse, is shown in Fig. 21a for the clamped beam covered by 0.015-in. 
thick HE layer (DuPont EL 506D) over a 2-in. span. The finite-span HE edge 
effect persists to a distance of about 0.5 to 1.0 inch. 

In theoretical analyses the impulsive loading can conveniently be ap- 
proximated by assuming sin initial velocity distribution. Corresponding to the 
present finite-element approach, the initial conditions to be specified are 
those nodal generalized initial velocities {q*} t=0 . From the spanwise experi- 
mental impulse distribution data shown in Fig. 21a, a uniform initial transverse 
nodal velocity is assumed to occur at those nodes of the beam elements which are 
entirely covered by the HE layer. However, for nodes within the HE edge-effect 
zone, the specification of initial velocities poses some uncertainty. Because 
the compatibility conditions required by the Bernoulli-Euler finite-element 
displacement model are that at boundary nodes of each element, the compati- 
bility of w and ij; (=&w/3r) - v/R) with neighboring elements is required, the 
initial velocities of node n (if it is located at the middle of the HE edge 
zone) may be specified by either (1) w = a, V = 0 or (2) w = a/2,V = 0, or (3) 
w = a/2, f = a/&, where "a" is the uniform initial normal velocity assigned to 
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the nodes covered by the HE layer but not in the HE edge zone of span "Z". 

These three initial velocity simulations of the effect of a finite span of the 
HE layer on the distribution of imparted impulse (as depicted in Fig. 21b and 
designated as IV1, IV2, and IV3) have been tested in the present analysis. 

The clamped-beam midspan deflection responses resulting from using these 
three different initial conditions are compared with each other and with experi- 
ment in Fig. 22. As shown in Fig. 22, it can be seen that the use of IV1 
initial velocity representation gives a response with higher amplitude than 
both that of experiment and those responses predicted by using XV2 and IV3, 
whereas the responses obtained by employing IV2 and IV3 initial velocity repre- 
sentations are close to each other and are in good agreement with experiment. 
Hence, it may be concluded that the IV1 initial velocity representation gives 
a higher amount of impulse than the actual impulse imparted to the beam, where- 
as the IV2 and IV3 initial condition representations tend to simulate the ex- 
perimentally-imparted impulse better than the first one does. In this example, 
10 uniform elements are used to model the half-span of the beam. But, when a 
coarser mesh is used, as shown in Fig. 23 for the free-ring central-plane 
separation responses, the IV2 and the IV3 initial condition simulations can 
give very different responses. However, both predictions approach each other 
when the mesh sizes become finer. From the above examples, it may be concluded 
that a finer space mesh would be required, especially near the edge of the HE 
layer, before a reasonably accurate representation of the initial impulse 
loading conditions can be obtained by the present computational method, as 
correlated with experiment. 

5. 3.6.2 Effects of the Humber of Spanwise and Depthwise 
Gaussian Integration Points 

Concerning the numerical evaluation of the integrals for determining 
the element properties {p} and [h] of Eq. 4.37, Gaussian quadrature has been 
employed to carry out the spanwise integrations over the length of the element 
and depthwise at each spanwise Gaussian point. Gaussian quadrature has been 
used also to evaluate the stress resultants (axial force, moment, and/or shear 
force) . 
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As illustrated by the clamped-beam example, the midspan-deflection 
responses obtained by using various numbers of spanwise Gaussian points are 
shown in Fig. 24; in each case, 4 depthwise Gaussian stations have been em- 
ployed to evaluate the stress resultants at each spanwise Gaussian point. It 
is seen that the prediction obtained by using only one spanwise Gaussian point 
tends to deviate appreciably from the behavior predicted by using 2, 3, or 4 
spanwise Gaussian points; the 2-spanwise-Gaussian-point result tends to be 
somewhat too stiff; while the 3- and 4-spanwise-Gaussian-point results are 
very close to each other. Also shown in Fig. 24 are the results obtained by 
assuming that the stress resultants may be approximated over the length of the 
element by their values at the center of the element; this prediction is seen 
to be better than that for the 1-spanwise-Gaussian-point case, but the 
"structure" tends to be too stiff. 

As for the effects of varying the number of depthwise Gaussian points to 
evaluate the stress resultants, shown in Fig. 25 are the midspan-deflection 
responses of the clamped-beam obtained by using 2, 3, 4, and 5 depthwise 
Gaussian points. In each case, 3 spanwise Gaussian points are used. It is 
seen that there is not very much difference among the deflection responses, as 
the number of depthwise Gaussian points is increased from 2 to 5. This is 
probably because the stretching behavior is predominant for the present clamped- 
beam example. The responses of the 2- and of the 4-point case differ somewhat, 
while the use of more than 4 depthwise Gaussian points affected the predicted 
response only very little. + In view of the above results and those of Refs. 44 
and 48, it appears reasonable to conclude that the use of 3 spanwise Gaussian 
points (or stations) and 4 depthwise Gaussian points at each spanwise Gaussian 
station suffices for (a) representing the internal stress distributions across 
the elements thickness and (b) the spanwise integration over the element length. 

5. 3.6. 3 The Effect of Strain- Rate Sensitivity 

In order to illustrate the effect on the transient response of using 
material strain-rate sensitivity in the present analysis, the impulsively- loaded 

+ Sirailar calculations for the free ring (Fig. 11c) using 3 spanwise Gaussian 
stations and either 4 or 6 depthwise Gaussian stations exhibited very little 
difference in the predicted responses. 
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clamped beam has been analyzed with the material property approximated either 
by EL-PP (elastic-perfectly plastic) or by EL-PP-SR (elastic perfectly-plastic 
and strain-rate sensitive) . The midspan deflection responses are compared in 
Fig. 26. The EL-PP-SR solution gives an 3% reduction in peak deflection com- 
pared with the EL-PP results and this peak occurs about 40 dsec earlier than 
in the EL-PP solution. 

For the impulsively- loaded free ring with CC Bernoulli-Euler-type ring 
elements , the predicted centerline midplane separation responses obtained by 
using the EL-SH-SR (elastic, linear-strain-hardening, and strain-rate sensitive) 
approximation for the material behavior leads to a 29% smaller peak amplitude 
and this peak occurs about 400 yu sec earlier than in the EL-SH solution, as 
shown in Fig. 27. Shown also are the responses predicted by using LC Bernoulli- 
Euler-type ring elements with or without including the strain-rate effect. 

Figure 23 shows the strain-rate effect on the central deflection re- 
sponses of the impulsively-loaded clamped ring which was modeled by using 
either the LC Bernoulli-Euler-type element or the Timoshenko- type element. It 
is seen that including tne strain-rate effect produces a "stiffer response"; 
i.e., small peak deformation response and earlier time-to-peak compared with 
the corresponding strain-rate independent predictions. 

The above-mentioned peak deformation response reductions and earlier 
peak responses caused by assuming the material to be strain-rate sensitive 
(and to follow Eq. 2.75) is also observed in the finite-difference calculations 
of Ref. 44. 

5. 3.6.4 The Effects of Including Transverse Shear Deformation 

In order to examine the influence on the predicted response by using 
the present Timoshenko- type element as developed in Subsection 4.4 (which 
takes the transverse shear deformation into account) as compared with that 
obtained by using the LC and CC Bernoulli-Euler elements , the impulsively- 
loaded ring problems were analyzed with only the linear assumed displacement 
functions for v , w, and 0 for the Timoshenko-type element. 

The maximum natural frequency (linear system) as computed by employing 
the Timoshenko-type elements is compared with that obtained by using 
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Bernoulli -Euler-type elements, in Fig. 29 and Fig. 30 for the clamped beam 
an<_ the free ring, respectively, as a function of the number of elements per 
half-span. It is observed for both the beam and the ring that the use of 
Timoshenko-type elements gives a larger maximum natural frequency , and hence 
a small critical At, than that obtained by the use of Bernoulli-Euler-type 
elements . 

For the clamped beam, shown in Fig. 31 are the Timoshenko- type predic- 
tion and the Bernoulli-Euler-type prediction of the midspan deflection re- 
sponses. Good agreement between these solutions is observed. 

Turning to the clamped ring results, the predicted central deflection 
responses for the Timoshenko- type element are compared with those obtained by 
using (1) LC Bernoulli-Euler-type elements and (2) CC Bernoulli-Euler-type 
elements in Fig. 32. The observations are that the agreement of the Timoshenko' 
type element prediction with the CC Bernoulli-Euler-type prediction is far 
better than with the LC Bernoulli-Euler-type element prediction. 

It should be noted that, for the present beam and ring examples, the 
transverse shear deformation effect essentially can be neglected, because of 
the thinness of the beam (thickness/span * 0.0125) and ring (thickness/radius «= 
0.042) . Hence, these examples permit confirming: (1) that the present 
Timoshenko- type element provides accurate large-deflection elastic-plastic 
transient response predictions of Bernoulli-Euler-type deformation and (2) the 
deficiency in the LC-type element, but do not provide a critical evaluation of 
the present Timoshenko-type element to predict large-deflection elastic-plastic 
responses with significant transverse shear deformation effects. An appro- 
priate such example having a reliable solution (or test result) has not been 
located . 

5.3,7 Comparison of Accuracy and Efficiency of Finite-Element 
Solutions Versus Finite-Pifference Solutions 

5. 3. 7.1 Scope of Comparisons 

In Ref. 44, experimental measurements of transient deformations and 
strains for impulsively- loaded clamped beams, clamped rings, and free rings 
which undergo large-deflection, elastic-plastic responses are compared with 
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finite-difference (FD) predictions — wherein finite differencing is employed 
for both spatially- and time-varying quantities. Good theoretical-experimental 
agreement has been demonstrated. 

Since the free-ring example (Fig. lie) embodies the most reliably de- 
fined boundary conditions of the above-cited three cases, only this example is 
used (see Subsection 5.3.7. 2) to compare the present finite-element (FE) pre- 
dictions with FD predictions, and with experiment. Similar comparisons, 
carried out for the other two examples of Fig. 11, show similar comparative 
results. 

This free-ring example is also used in Subsection 5. 3.7.3 to illustrate 
and assess the comparative efficiency of the FD and the present FE prediction 
method (s) — in terms of the amount of computer central processing unit (CPU) 
time required to carry out calculations for a given time of actual structural 
response and at the same time to provide peak deformation (and/or peak strain) 
predictions within a given percentage of the converged value (displacement or 
strain) . This comparison is believed to provide a reasonably good , although 
tentative, assessment of the comparative cost for providing predictions of a 
"given accuracy" by the FE and the FD approach. 

5. 3. 7. 2 Comparison of Experiment with FE and FD Predictions 

The geometry, material properties , and loading conditions (represented 
here as initial velocity conditions) are shown, in Fig. 11c for the free ring. 

It has been demonstrated both in Ref. 44 and in Subsection 5. 3.6. 3 that the 
neglect of strain-rate effects in representing the mechanical properties of 
this 6061-T6 aluminum alloy material leads to a vast overprediction of the 
structural response. Accordingly, for convenience, the only material property 
representation employed in the present comparison is EL-SH-SR (elastic, linear 
strain hardening, and strain-rate dependent), with D • 6500 sec 1 and p “ 4 
(see Eq. 4.50). Also, the free ring is assumed to undergo Bernoulli-Euler-type 
deformation. 

In both the FE and the FD predictions being discussed, the temporal 
3-point central-difference operator is used. In view of the FE results dis- 
cussed in the previous subsections, only the improved formulation type of FE 
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predictions is included here; also, the following features are used as being 
appropriate and "adequate": 

(a) in each element, three spanwise Gaussian stations are used 
to evaluate {p} and [h). 

(b) four depthwise Gaussian stations are used to evaluate 
the inplane stress resultants and the moment resultants 
at each spanwise Gaussian station. 

(c) Bernoulli-Euler finite-elements having the CC type of 
assumed displacement function are employed. 

For the FD predictions, the method of Ref. 44 as subsequently improved 
and embodied in the JET 2 computer program of Ref. 169 has been used. In this 
method, the stress and moment resultants are evaluated only at each "finite- 
difference mass-point" station; at each such station, four depthwise Gaussian 
stations are used to carry out these evaluations in order to provide appropri- 
ate correspondence to the evaluations used in the present FE calculations. It 
should be noted that in the finite-difference calculations for the present 
type of structure, there are only two degrees of freedom (axial displacement 
and lateral displacement) at each space-mesh intersection (also called "mass 
point station"). Also, the mass matrix in the FD method is obtained by 
(automatic) lumping. 

If one were to take advantage of symmetry, one could model the half ring 

by a number N of finite elements of the CC- type assumed-displacement function 
E 

or by N q finite-difference space-mesh stations. Accordingly, the associated 
number of degrees of freedom for each would be as follows: 


Number 

Degrees of Freedom 

Unrestrained (dof ) 
u 

Retained Number of 
Degrees of Freedom 
with Symmetry Re- 
straints Applied (dof^). 

FE: N 

4N_ + 4 

4N_ 

E 

E 

(for CC elements) 

E 

FD: N 

2N_ 

2N_ 

D 

D 

D 


In order to illustrate typical comparisons of the FE vs FD predictions. 
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calculations which utilize a roughly comparable number of degrees of freedom 

are selected for presentation here. Accordingly, predictions for N* = 18 (or 

+ “ 
dof^ = 76 and dof g = 72) and N Q = 40 (or dof^ = dof g = 80) are compared in 

Fig. 33a for the centerline midplane separation history and in Fig. 33b for 
the deformation profiles at 1140 usee and 2580 usee. It is seen that there is 
very good agreement between the FE and the FD prediction shown here; also, both 
predictions are in reasonably good agreement with experiment. Note that the 
experimental deformed-ring profiles exhibit some asymmetry, possibly from 
initial out-of-roundness and/or some unintentional and undefined asymmetry in 
the applied impulsive loading. Whereas, in both the FE and FD predictions, 
symmetry was imposed by choice — nonsymmetric cases could be analyzed, how- 
ever. 

For the improved- formulation FE analysis of this ring which was modeled 

by 18 CC elements, transient response calculations were carried out by using 

(a) the consistent mass matrix and (b) the lumped mass matrix. In each case, 

the At used was approximately 0.8 At^ * 0.8 [2/( d^ ay ] as numerical experiments 

had previously verified to be acceptable; accordingly, the At's were, 0.6 and 

1.8 usee, respectively. Since to the scale of the present plots these two 

transient response predictions were nearly alike, only the "more efficient" 

lumped-mass FE prediction is shown in Fig. 33. For the FD calculation, the At 

1/2 

employed was 9/8 usee which is 99% of ^critical = (As) /(E/p) where As is 
the finite-difference mesh length. 

An examination of a more sensitive quantity is provided by the dynamic 
strain responses. Predicted and experimentally-measured strain-time histories 
at 6 locations on the ring are compared in Figs. 33c and 33d. Again, both the 
FE and the FD predictions are in reasonably good amplitude and phase agreement 
with experiment, with the FE results being somewhat better. 

Strain profiles predicted and measured at the 3,000 usee instant are 
compared in Fig. 33e. Fairly good agreement among the FE prediction, the FD 
prediction, and experiment is observed. Also, note in Fig. 33e that abrupt 
reversal of the strain occurs at 8 = 60° which is the location of the edge of 
the high explosive layer, and that the ring undergoes essentially pure bending 
for 8 greater than about 70°. Considerable compression strain plus bending 

+ Uniform element lengths and uniform mesh lengths were used. 
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strain is seen to exist for small-6 locations. 


5. 3. 7. 3 Convergence and Efficiency Comparisons 

It was noted in Subsection 5.3.3 that the conventional finite-element 
formulation requires more computer storage and running time to analyze a given 
problem than does the improved finite-element formulation. For example, for 
1000 cycles of computing for the clamped-beam example of Fig. 11a, the con- 
ventional formulation required 1.85 minutes while the improved formulation re- 
quired 1.41 minutes of CPU time (a saving of about 24%) where the same number 
of LC elements was used for both computations; the predictions from these two 
calculations were almost indistinguishable. Since it is clear that the improved 
formulation provides more efficient predictions , only the improved formulation 
is used in subsequent comparisons in this subsection. 

For assessing the comparative efficiencies of the FE approaches (a) con- 
sistent mass and (b) lumped mass versus the finite-difference calculation, the 
free-ring example of Fig. 11c is used for both convenience and the fact that 
large deformations and elastic-plastic transient responses are involved. 

By increasing the number of equal-length finite elements (of the CC type) 
to model the half ring (taking advantage of symmetry) , transient response pre- 
dictions have been carried out. The following two useful indices of the re- 
sponse were monitored, and are discussed herein, in obtaining a useful measure 
of convergence; 

(1) the peak relative displacement of the ring at the symmetry 
plane (or (£_ ) . 

(2) the peak circumferential strain at several locations on 
the outer surface of the ring. 

2 

By plotting these results as a function of (1/N ) , and extrapolating to 

2 ■ 

(1/N ) = 0, "converged" results were estimated. The peak responses predicted 

were then ratioed to the appropriate "converged" result. Accordingly, shown 

in Fig. 34a is the ratio to the converged result of the predicted peak center- 

line relative displacement as a function of the number of unrestrained degrees 

of freedom dof (= 4N + 4). Shown in Fig. 34b is the average at 6 = 87°, 

U E 
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92.5°, and 177® of the peak circumferential strain ratioed to the respective 
converged value as a function of dof^; a mean dashed line is shown as well as 
a speckled band to indicate that such results show scatter depending upon the 
8-location chosen and the number of 8-locations which one could use to con- 
struct this "average". It should be noted that, as discussed in Subsection 
5. 3. 6.1, the specification of initial nodal velocity at the node of the finite 
element which is located at the spanwise edge of the high-explosive layer poses 
some uncertainty; that is, the use of different initial nodal velocities 
(IV2 or IV3) will yield very different coarse-mesh deformation responses be- 
cause of the nature of the initial-velocity distribution. Improvements earn be 
made by using finer meshes near the edge of the high-explosive layer. However, 
pending a more rational way of representing the finite span impulsive loading 
can be devised, the results presented in Figs. 34a and 34b are based on a 
uniform-mesh size and the IV3 type of initial-velocity representation. 

Similar calculations and reductions were carried out for the finite 
difference method using the JET2 computer program. Shown in Fig. 35a is the 
ratio to the converged value of the predicted peak centerline relative dis- 
placement as a function of the number of degrees of freedom (i.e., 2N q ) . 
Similarly, Fig. 35b shows the average at 8 = 27°, 45°, 63°, 81°, 99° and 171° 
of the peak circumferential strain ratioed to the respective converged value 
as a function of the number of degrees of freedom in the finite-difference 
model; for these 8-locations, the peak strain ranged from about 2 to 4 percent. 
In Fig. 35b a dashed "average" curve through a speckled band is shown to indi- 
cate that such results show scatter which depends upon the 8-location studied. 

With respect to computing time, the amount of computer central process- 
ing unit time is believed to provide a reasonable basis for comparing the 
"time consumption" of the FE and the FD calculations. For the improved formula- 
tion version of the finite-element method when applied to large-deformation 
elastic-plastic transient response problems, it has been found that the maximum 

allowable time-increment size At is about 0.8(2/u } where u) is the largest 

max max 

natural frequency contained in the finite-element model of the structure for 
small vibrations, where the 3-point central-difference finite-difference time 
operator is used. When the present CC type of finite element is used together 
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with (a) the consistent mass (C) matrix or (b) the lumped mass (L) matrix, the 

u> is found to be (see Fig. 16a) : 
max * 

(u>_) = [-0.1445 + 0.03578(dof )]10 6 radians/sec (5.8a) 

c max u 

( 0 ) ) = [-0.0927 + 0.0127 (dof )]10 6 radians/sec (5.8b) 

L max u 

where dof = 4N + 4 for the half ring; note that (gj ) « (w„) 

u E 9 L max C max 

Further, it has been found that the amount of central processing unit time 

(CPUT) in minutes , for the improved FE method, is given approximately by: 

CPUT = [42 x 10 _6 ] [dofj (minutes) (5.9) 

where t is the number of seconds of actual structural response to be computed 

and At is also in seconds. Since the allowable At * 0.8 (2/u ) , it follows 

max 

that: 

(CPUT) C = [-3.794(dof u ) +0.9393 (dof u ) 2 ]t (5.10a) 

(CPUT) = [-2.433 (dof ) +0. 3325(dof ) 2 ]t (5.10b) 

X* u u 

Similarly for the JET 2 program which uses the finite-difference method, 
it has been found that 

(CPUT) fd = [23 x 10" 6 ] [dofj (minutes) (5.11) 

where dof^ = 2N^. Also since the largest allowable At is given by 

(As)/(E/p)^ 2 , it follows for the half ring with R = 3 in. and At for conserva- 

1/2 

tism taken as 0.99 (As) /(E/p) that 

(CPUT) = 0.2513 (dof ) 2 t (5.12) 

FD U 

With the "convergence" results of Figs. 34a through 35b and with the 
above central processing unit time for computing , one may estimate the compara- 
tive CPUT values for the FE_ , FE , and FD computer programs used here in order 
to predict the peak centerline relative deformation or the peak circumferential 
strain to within selected percentages of the converged value for each type of 
calculation. An example of such estimates is tabulated as follows for 
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computing 1500 usee of transient response of the free ring, wherein it has 
been assumed that the "convergence rate” for the lumped mass version of the 
finite-element method is essentially the sane as that shown in Figs. 34a and 
34b for the consistent mass version of the finite-element method. 


CENTRAL PROCESSING UNIT COMPUTING TIME ON IBM 370/155 
AT MIT (MINUTES) FOR 1500 MICROSECONDS OF ACTUAL 
STRUCTURAL RESPONSE 


Finite-Element 
Improved Formulation 

Finite Difference 


CC and 

Consistent Mass 

CC and 
Lumped Mass 


Percent of "Converged" 
Peak Ring Relative 

Displacement 




3 

2.03 

.65 

1.36 

2 

3.24 

1.06 

1.64 

1 

6.89 

2.32 

2.01 

Percent of "Converged" 
Peak Strain on the 
Average 




10 

2.48 

.82 

1.85 

5 

5.58 

1.87 

2.47 

3 

7.50 

2.53 

2.92 


125 
















Comparing the two types of finite element calculations (consistent 
mass and lumped mass) , it is evident that it is substantially more efficient 
to use the lumped-mass version of the finite element method. Further, since 
the transient responses of the consistent mass and of the lumped mass FE cal- 
culations differ only slightly, it is recommended that the lunped mass version 
of the FE method be selected as being more efficient and adequate for engineer- 
ing prediction purposes. 

Note also that this example comparison indicates that the FE. calcula- 

ii 

tion is often more efficient than the FD prediction method. However, the con- 
clusion that the FE^ calculation will usually be more efficient than the FD 
method is not warranted on the basis of this limited comparison. Many more 
examples and much more thorough comparisons would be required before any judg- 
ment of this type could be made — it would not be unexpected to find that one 
method would be superior for certain types of examples and the other method 
would be superior for certain other types of structural transient response 
problems . 

Finally, in comparing FE T calculations with FD calculations wherein each 

id 

employs the same number of degrees of freedom dof^, for the free ring, it is 
seen that the ratio (CPUT) __/ (CPUT) t for dof = 10, 50, and 100 is 2.815, 0.884, 
and 0.816, respectively, and asymptotically approaches 0.755 as dof^ is in- 
creased indefinitely if it is presumed that Eqs. 5.10b and 5.12 would apply — 

these equations would, of course, cease to be valid where dof becomes so 

/ u 

large that resort to auxiliary storage and retrieval would be needed. The 
reasons for this relative computing time consumption involve the facts that: 

(a) For the same number of degrees of freedom, the total number 
of spanwise Gaussian points employed in the FE method is 1.5 
times as large as the total number of space-mesh intersections 
used in the FD method; Gaussian evaluations in the latter are 
performed only at the space-mesh intersection stations. At 
each spanwise Gaussian station, inplane stress resultants and 
moment resultants are evaluated. 
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(b) In the FD method, the strain is computed from the strain- 
displacement relations by finite-difference approximations 
in terms of displacements at two or three neighboring space- 
mesh stations. On the other hand, in the FE evaluation of 
strains, one computes the strains from the strain-displace- 
ment relations directly by using the assumed-displacement 
field of the finite element involved; this involves a time- 
consuming matrix multiplication. Of course, one could cir- 
cumvent this in part in the FE method by resorting to the 
use of finite-difference approximations for each spanwise 
Gaussian station in terms of nearby nodal generalized dis- 
placements . 

Also, several recent papers (Refs. 187, 138, and 189) in which various 
aspects of the finite-difference method versus the finite-element method are 
discussed have just appeared; these documents are recommended reading — cover- 
ing some of the present aspects as well as others. It perhaps should be men- 
tioned that the finite-difference equations formulated in Refs. 39, 188, and 
189 are based upon the variational-energy principle; the derivatives of the 
field variables in the variational functional are replaced by appropriate 
finite-difference quotients which involve only the values of the variables at 
the space-mesh stations. This finite-difference formulation is, hence, somewhat 
different from that described in Section 3, wherein the derivatives in the 
governing differential equations and boundary conditions are replaced by appro- 
priate finite-difference expressions. 
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SECTION 6 


APPLICATION: CONTAINMENT/DEFLECTION RING RESPONSES 
TO ENGINE ROTOR BLADE FRAGMENT IMPACT 

6.1 Introduction and Problem Definition 

Since the advent of the turbojet engine, there have been, from time to 
time, failures of turbine and/or compressor rotor blades and/or disks on en- 
gines of both military and civilian aircraft (Refs. 1, 170-173). Fragments 
which are uncontained (that is, penetrate the engine casing) might injure 
personnel occupying the aircraft and might cause additional damage to fuel 
lines and tanks, control systems, and other vital components . Although strenu- 
ous efforts have been and continue to be made to avoid blade/disk failures 
through improved materials, design, fabrication, and inspection, a not-insignifi 
cant number of such failures persist. It is desirable, therefore, to provide 
protection (a) for on-board personnel of aircraft in flight and (b) for vital 
components . 

Similar but perhaps less severe fragment containment/control problems 
may be encountered where turbines and/or energy-storing flywheels are used in 
stationary power plants, aboard ships, and/or in land vehicles such as buses, 
trucks, automobiles, and train locomotives. In these cases, there would usually 
be less concern about the "weight penalty" for insuring fragment containment 
than for aircraft. 

Two distinct avenues for providing this protection are evident. First, 
the structure surrounding the "failure-prone" rotor region could be designed 
to contain (that is, prevent the escape of) the rotor-burst fragments completely 
Second, the structure surrounding this rotor could be designed so as to prevent 
fragment penetration in, and to deflect fragments away from, certain critical 
regions or directions but to permit fragment escape readily in other "harmless" 
regions or directions. One or both of these schemes could, in principle, be 
employed in a given design. In any event, this desired protection is sought 
for the least weight and/or cost penalty. If only one of these two schemes 
were to be adopted, one might expect that the second would be most cost and/or 
weight effective. However, the present (1) knowledge of the fragment-structure 
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interaction phenomena and (2) analysis/design tools are inadequate to permit 
making a definitive comparative assessment at this time, although much progress 
has been made in both of these areas in the past few years. 

As pointed out in Ref. 1, NASA has been sponsoring a research program 
which is designed to meet the objective of providing the necessary protection 
to aircraft without imposing large weight penalties. Starting about 1964, the 
Naval Air Propulsion Test Center (NAPTC) under NASA sponsorship has constructed 
and employed a spin-chamber test facility wherein rotors of various sizes can be 
operated at high rptn, failed, and very importantly the interactions of the re- 
sulting fragments with various types of containment and/or deflection structures 
caui be studied with high-speed photography, in addition to post-mortem studies 
of the containment/deflection structure and the fragments. Many such tests in- 
volving single fragments or many complex fragments impinging upon containment 
structures of various types and materials have been conducted (Refs. 170-174) 
and have substantially increased the body of knowledge of the attendant phe- 
nomena. For the past several years NASA has sponsored a research effort at 
the MIT Aeroelastic and Structures Research Laboratory (ASRL) to develop methods 
for predicting theoretically the interaction behavior between fragments and con- 
tainment/deflection structures, as well as the transient deformations and re- 
sponses of containment/deflection structures — the principal objective being 
to devise reliable prediction/design procedures and containment/deflection 
techniques. Important cross-fertilization has occurred between the NAPTC ex- 
perimental and the MIT-ASRL theoretical studies, with special supportive- 
diagnostic experiments and detailed measurements being designed jointly by 
NASA, NAPTC, and MIT personnel and conducted at the NAPTC. Subsequent analysis 
and theoretical-experimental correlation work has been increasing both the under- 
standing of the phenomena involved and the ability to predict these interaction/ 
structural-response phenomena quantitatively. 

Because of the multiple complexities involved in the very general case 
wherein the failure of one blade leads to impact against the engine casing, re- 
bound, interaction with other blades and subsequent cascading rotor-failures 
and multiple-impact interactions of the various fragments with the casing and 
with each other, it is necessary to focus attention initially upon a much 
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simpler situation in order to develop an adequate understanding of these col- 
lision-interaction processes. Accordingly, rather than considering the gen- 
eral three-dimensional large deformations of actual engine casings under multi- 
ple rotor-fragment attack, the simpler problem of planar structural response 
of containment structures has been scrutinized. That is, the containment 
structure is regarded simply as a structural ring lying in a plane; the ring 
may undergo large deformations but these deformations are confined essentially 
to that plane. For such a case, a numerical method of analysis to predict the 
transient large-deformation responses of such structures to known impulsive 
and/or transient external loading has been developed at the MIT-ASKL and has 
been verified (Refs. 44 and 45) by evaluative comparison with high-quality ex- 
perimental data, to provide reliable predictions. This prediction method is 
sufficiently simple that one can feasibly carry out certain types of parametric 
structural response calculations, provided that known or prescribed externally- 
applied forces or impulses are employed; limited such studies are reported in 
Refs. 175 and 176. 

In the present context, therefore, the crucial information which needs 
to be determined (if the structural response of a containment ring is to be 
predicted reliably) concerns the magnitude, distribution, and time history of 
the loading which the ring experiences because of fragment impact and inter- 
action with the ring. Two means for supplying this information have been 
considered : 

(1) The TEJ concept (Refs. 169, 176, 177) which utilizes measured 
experimental ring position-time data during the ring-fragment 
interaction process in order to deduce the external forces 
experienced by the ring. This concept has been pursued. An 
important merit of this approach is that it can be applied 
with equal facility to ring problems involving simple single 
fragments such as one blade, or to cases involving a conylex 
multi-bladed-disk fragment. The central idea here is that if 
the TEJ type analysis were applied to typical cases of, for 
example, (a) single-blade impact, (b) disk-segment impact, 
and/or (c) multi-bladed disk fragment impact (see Figs. 36 
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and 37) , one could determine the distribution and time 
history of the forces allied to the containment ring 
for each case. Such forces could then be applied tenta- 
tively in computer code response-prediction-and-screen- 
ing studies for similar types of ring-fragment interaction 
problems involving various other materials, where guidance 
in the proper application of these forces or their modifica- 
tion could be furnished by dimensional-analysis considera- 
tions and selected spot-check experiments. 

On the other hand, this approach suffers from the fact that ex- 
perimental transient structural deformation data must be avail- 
able; the forcing function is not determined from basic material 
property, geometry, and initial impact information. 

(2) The second approach, however, utilizes basic material 
property, geometry, and initial impact information in 
an approximate analysis which employs the basic prin- 
ciples of energy and momentum conservation as well as 
material property constitutive data. If the problem 
involves only a single fragment, this method can be 
carried out and implemented without undue difficulty, 
but can become very complicated and time consuming if 
complex fragments and/or multiple fragments must be 
taken into account. 

Approach 1 is explained in detail in Refs. 169, 176, and 177. The present 
report deals with approach 2 and confines attention to problems involving only 
a single simple fragment; problems involving more complicated -fragments are 
left for future consideration. 

Various levels of sophistication may be employed in approach 2. One 
could, for example , employ finite-difference methods wherein both the contain- 
ment ring and the fragment are represented by a suitably fine three-dimensional 
spatial mesh and the conservation equations are solved in time for simple con- 
figurations by digital computer codes such as HEMP (Ref. 178), STRIDE (Ref. 179) , 
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and/or HELP (Ref. 180)/ which take into account elastic, plastic, strain harden- 
ing, and strain-rate behavior of the material. Such computations , while vital 
for certain types of problems, are very lengthy and expensive, and are not well 
suited for the type of engineering analysis/design purposes needed in the present 
problem; for complicated or multiple fragments, such calculations would be pro- 
hibitively complicated, lengthy, and expensive. A simpler, less complicated , 
engineering-analysis attack within this general framework is needed. 

Two categories of such an engineering analysis may be identified and 
are termed: (a) the collision- imparted velocity method (CIVM) and (b) the 
collision- force method (CFM) . The essence of each method follows: 

(a) Collision- Imparted Velocity Method (C1VM) 

In this approach the local deformations of the fragment or of the 
ring at the collision interface do not enter explicitly, but the 
containment ring can deform in an elastic-plastic fashion by mem- 
brane and bending action as a result of having imparted to it a 
collision-induced velocity at the contact region via (a) per- 
fectly-elastic or (b) perfectly-inelastic behavior. In fact, 
any type of material behavior may be accommodated readily. 

Since the collision analysis provides only collision imparted 
velocity information for the ring and the fragment ( not the 
collision-induced interaction forces themselves) , this procedure 
is called the collision-imparted velocity method. 

(b) Collision-Force Method (CFM) 

In this method the primary information predicted in the collision 
analysis consists of the collision-induced interaction forces 
themselves ; the associated and subsequent ring and fragment re- 
sponses are also predicted. 

Since the CIVM is much simpler to implement than the CFM approach, the CIVM 
scheme has been studied and is discussed (later) in some detail in Subsections 
6.2, 6.3, and 6.4 of this report. The CFM approach, described briefly in 
Subsection 6.2.2, is currently under study, development, and feasibility 
evaluation; these findings will appear in a future report. 
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In Subsection 6.2 approximate collision analyses cure discussed in de- 
tail. Some applications and evaluation of the CXVM approach are documented 
in Subsection 6.3. Finally, some comments on potential extensions of the CIVM 
scheme to more complex problems are given in Subsection 6.4. 

6.2 Approximate Collision Analyses 

Under consideration in this subsection sure approximate methods for pre- 
dicting the "immediate consequences" of the impact of a fragment against 
another physical body. One may regard the fragment as being (a) rigid, (b) 
perfectly elastic, (c) perfectly inelastic, or (d) deformable elastic-plastic 
with EL-SH-SR behavior; similar behavior may be attributed to the body which 
is struck by the fragment. When both bodies are treated as behaving according 
to (d) , the modeled constitutive behavior of each most closely simulates the 
true physical behavior, but the associated impact-interaction is the most com- 
plex of the various options. For engineering applications purposes, one de- 
sires to employ the simplest and least expensive procedure which will give ade- 
quate engineering accuracy. Accordingly, various convenient and plausible 
assumptions are invoked to predict the "immediate impact- interaction" behavior 
of a fragment which collides with a containment/deflection structure. 

For present purposes, the subject collision- interaction problem is simpli- 
fied by restricting the motion to lie in one plane; the extension of the analy- 
sis to the more general three-dimensional motion-and-deformation behavior can 
be carried out, if desired, in a future investigation. Also, only a single 
simple fragment is considered, as depicted schematically in Fig. 36a. As 
noted earlier, principal attention is given in this report to the approximate 
analysis which is termed the collision-imparted velocity method (CIVM) ; only 
limited discussion is devoted to the collision-force method (CFM) in 
Subsection 6.2.2. 

6.2.1 Collision- Imparted Velocity Method (CIVM) 

For the CIVM approach, the following additional simplifying assumptions 
are invoked: 

(1) In an overall sense, the fragment is treated as being rigid. 

It does not undergo bending or extensional deformation, but 


133 



at the "immediate contact region" between the fragment 
and the struck object (termed the "target" for conveni- 
ence) , the collision process is regarded as being in- 
stantaneous with a perfectly-elastic or a perfectly- 
inelastic interaction. 

(2) The colliding surfaces of both the fragment and the target 
are perfectly smooth; hence, no forces and/or velocities 
(or momentum) are either transmitted or imparted in the 
tangential direction. 

(3) During the collision, the contact forces are the only ones 
considered to act on the fragment and (in an antiparallel 
fashion) on the ring. The internal forces are approximated 
as being zero because the duration of the impact is so 
short as to preclude their "effective development". 

(4) The collision process is instantaneous and involves only 
the fragment and the containment-ring segment which en- 
compasses the ring-fragment collision point, as indi- 
cated schematically in Fig. 38. The word "instantaneous" 
implies that the internal forces can be neglected be- 
cause stress waves have not had time to traverse the 
area near the collision from regions reasonably remote 
from the impact point in the ring. It is this in- 
stantaneity which permits one to omit the internal 
forces from the CIVM model. 

(5) To avoid unduly complicating the analysis and because 
of the smallness of the arc length of the target-ring 
element, the ring element is treated as a straight 
beam in the derivation of the impact equations. 

Two different approximation models are evident for the collision-interaction 
calculation for the fragment with the ring; namely, the consistent mass model 
and the lumped mass model. The former scheme treats each ring element as 
having a distributed mass; the latter scheme, employed in Ref. 181, considers 
the mass of each ring element to be concentrated at its two end nodes. 
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Accordingly, there are two sets p£ impact equations corresponding to these two 
approximate interaction models, which are described, respectively, in Subsec- 
tions 6. 2. 1.1 and 6. 2. 1.2. 


6. 2. 1.1 Consistent-Mass Collision Model 

In this collision model, the affected ring segment is idealized as a straight 
beam of length s having a nonuniform distribution of mass. Since perfectly 
smooth + surfaces for both the ring segment and the fragment are assumed to 
exist at the impact location, the instantaneous collision process results in 
an equal and opposite impulse applied to the ring (beam) segment and to the 
fragment in the direction normal to the axis of the ring segment; accordingly, 
tangential-component velocities are unaffected — only the normal-direction 
components of the velocities of these two bodies undergo- change because of 
the collision. Hence, only these components are employed in the following con- 
servation relations. 


Referring to the schematic and notation of Fig. 39a and to the idealized 
"line" geometry depicted in Fig. 40a, the impulse-momentum law and the kinetic 
energy conservation law may be written to characterize the "instantaneous im- 
pact behavior" of this system, as follows: 


Translational Impulse-Momentum Law 

m r ( ( U, - Ur ) + a < u)/ - “Jr)] = "f>„ 

7n f f U f '~ U f ] = - f„ 

Rotational Impulse-Momentum Law 

I r? [ - ^r] --f n (d-Y)S 

If [ (Of - u)f]= ( i Sin B ) 


(ring segment) 


(fragment) 


(ring segment) 


(fragment) 


( 6 . 1 ) 

( 6 . 2 ) 

(6.3) 

(6.4) 


+ 

The smooth surface assumption is invoked, for convenience, in the present 
analysis, but could be relaxed, if desired, in future work. 
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Conservation of Kinetic Ener gy 

2 , .1 


YY\y 


r [ U r A C0 r J i- j l rc , [ J + frlf [ ] + £- If [ 1 


- ^ >n r [ U r f a C0 r ] ■+ J Irf [ dr] 


+ T- + T Ii 


? J f ^ 


(6.5) 


where 


= mass of the ring element which has a length s 

m f = mass of the fragment whose length is + Jt f2 (Fig. 39a) 

SL in Eq. 6.4 = fragment length from the impact point to the c.g. 
of the fragment 

I = mass moment of inertia of the ring element 
about its center of gravity (c.g.) 


7YS , r'bS 2 

* r,f § ) § c/ § + ( m J ^ d 1 


(see Fig. 39a) 


m(£) or m(ri) denotes mass per unit length of the ring element 

I f = mass moment of inertia of the fragment about its c.g. 

p^ = normal impulse 

.1 

a = (- - Y )s 
U l + U 2 


° 1- U 2 


= angular velocity of the ring element 


U^,U 2 = normal velocities at ring-element nodes 1 and 2, 
respectively, immediately before impact 

, 0 )^ = fragment c.g. normal-direction velocity and angular 
velocity immediately before impact 
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u;, o)^, V. O*, U', <j)' £ = translational normal-direction velocities 

and angular velocities immediately after 
impact (primes indicate these after-impact 
quantities) . 


Equations 6.1 through 6.5 represent five equations expressed in terms of the 

five unknowns U', <u' , U', 0 )', and p . 

r r r r n 

Next, it is convenient to eliminate p^ by using the Eq. 6.1-value for 
P n to replace in Eqs. 6.2, 6.3, and 6.4 to obtain, after dividing each 
equation by ra^ : 

Uf “ U , = - yy\ [ ( ■+ a odr) - ( (J r + a )] ( 6 . 2 a) 


I, [ U)/ - Ldr ] - -y m [ ( 1//+ t lair )~ + ^ ) ] 


(6.3a) 


I 2 ( ~ rtf J - z m { ( U' ■* a. oV )~ ( U r ■+ £ oV ;J 


where 


m = 


m r 



y = (oi -Y ) 5 


Z = i 5/n 0 


i£ 

"’f 


(6.4a) 


( 6 . 6 ) 


Rewriting and dividing Eq. 6.5 (the kinetic energy equation) by ra f one obtains, 
with the use of Eq. 6.6: 

j m[ ( U/+ a. (A? ) Z - ( U K + a u) r f] + j- i ( [( u)//- iu) r ; Z J 

+ j r ( u f / ' )Z - + j- i z [ ( f- f / j = o 

1 (6.5a) 

If to the second, third, and fourth terms of Eq. 6.5a, one applies, respectively, 
Eqs. 6.3a, 6.2a, and 6.4a, one obtains: 
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J m j ( Ur ■+ a. U),f )- ( Ur + a ) Z j 
4 J I, { u)/+ U3 r ] [ = ^- { ( Ur'+a - ( (J, •+ a i*V )J j 
■+2' ( Uf -*■ Uf | — ft) ( ( Ur + A M/ ) ~ ( U r + & M r ) J } 

+ { h |“V 4£ ^} {^f [ ( u;+au>/)- (Ur + a^r))] = o 

lz J (6.7) 

Next, dividing Eq. 6.7 by y( (U^ + auT ) - (U^ + aw^) ] , 
one obtains 

j [ U r 'Mj-d)S <J- [ U/- 2 j=-J [ Ll r +^-«*)5ttVj-fU f -2^jj (6.8) 

Equation 6.8 states that the velocity at the ring-impact point (C) relative 
to the velocity at the fragment impact point (A) is singly reversed by per- 
fectly-elastic impact, since one may readily verify, for example, that 

U„ - (Jf -( 1 Sin 6 ) u)^ (6.9a) 

U c = (J r + ( y =-U l ~ h (U i ~U,)0l (6.9b) 

Thus, if desired, one may express Eq. 6.8 by the self-evident notation 

= -\u c -u A ) 

As pointed out in Ref. 182 (pg. 4-37) , experiments on direct central 
impact of spherical bodies have shown that the relative velocities of spheres 
after impact are always less than before impact , and that these relative ve- 
locities are opposite in direction. The ratio of the relative velocity after 
impact to that before impact is called the "coefficient of restitution” and 
is generally denoted as e. It is found that 0 <_ e <_ 1, where e * 1 represents 
a perfectly-elastic impact and e *= 0 denotes a perfectly-inelastic impact. 
Typically, for glass e is 15/16, for ivory 8/9, steel and cork 5/9, wood about 
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1/2, and clay and putty 0.* Further, it should be recalled that during imper- 
fect impact of two colliding bodies of masses m^ and m^, with initial velocities 
and V^» there is a kinetic energy loss given by 


K.E.Loss 


m, m i 
l ( rn,+ wi 2 ) 


( V, -Vz )* U-e ) 


(6.10) 


For the perfectly-elastic case of e * 1, no kinetic energy loss occurs. For 
e / 1, there is a loss of kinetic energy; however, since the total energy 
must be conserved, the "lost kinetic energy" is simply converted to other forms 
such as thermal or heat energy, etc. A proper accounting of all of the energy 
could be done in a complete thermo-mechanical analysis. For the present ap- 
proximate analysis, however, one need not keep account of this kinetic energy 
loss. 


Applying the concept of the coefficient of restitution, e, to Eq. 6.8, 
one may "generalize" this equation to read: 

- r U'- H } = - e j [ U t +(a-y) -[U f -2 «*] } (e.iu 

where (1/2 - ot) s = (a - y) has been used. 


Finally, one can solve Eqs. 6.2a, 6.3a, 6.4a, and 6.11 to obtain ex- 
pressions for LP , , U|, and u)^, the "unknown" after- impact quantities. First, 

solving Eq. 6.3a for in terms of 


U) 


/ 

Y 


c0 t 


y m 

1,+yam 


( U,' - (J r ) 


(6.3b) 


Next, 


applying Eq, 6.3b to Eq. 6. 2a and solving for in terms of IT : 

u; = U f - ( US- Ur) 


(6.2b) 


Also, applying Eq. 6.3b to Eq. 6.4a and solving for in terms of U^: 
/ I, m Z . / . 

+ 1 (1,+yantj ( Ur~U r ) 


(6.4b) 


It will be shown subsequently that predicted containment ring structural re- 
sponses to fragment impact are insensitive to the value of e employed for 
0 < e < 1. 


139 



Finally, applying Egs. 6.2b, 6.3b, and 6.4b to Eg. 6.11, the following ex- 
pression for U* is obtained: 



U r 


J_+ 

K 



(6.12) 


Thus, using Eg. 6.12, Egs, 6.2b, 6.3b, and 6.4b yield., respectively. 


UL = 


Ur + U R 

t I t + yaw l\c * 


( 6 . 13 ) 


LOyf = WJf. + 


y m 


I'+yam K c 


) u. 


(6.14) 


0)r = U) 


I,wE 




f J tl+yam) K c 


(6.15) 


Since U* * (U ’ + U‘ )/2 and u>‘ * (O' - U*)/s, and similarly for the unprimed 
i a * r x * 

preimpact guantities, it follows that 



, y ” 

S 

l / 

i+e 

') Ur 


u,-{ 

1 1 , + yctw 

z 

J ( 

Kc 

( 6 . 16 ) 


1 -4 ym 

s 

K 

|+€ 

) u R 


U z j 

I ( +ya»i 

2 

Kc 

( 6 . 17 ) 


where 

K c * I + I( 7y7m f i; 2 1 (6 - 18a> 

u r' f U r + - [ ll f - Z 

= r p u, u z ] -t (6.18b) 


= U c - U A = relative velocity of the impact points: C and A. 

Note that subscript C of K c denotes that this guantity is associated with the 
consistent-mass collision model. 
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6. 2. 1.2 Lumped-Mass Collision Model 

For this collision calculation, the ring segment is treated as having 
only point masses and m 2 at nodes 1 and 2, respectively, rather than a 
distributed mass, as indicated schematically in Figs. 39b and 40b. The im- 
pulse-momentum law and the kinetic energy conservation law can be applied to 
this model also to obtain: 


Translational Impulse-Momentum Law 

(w, + ffl i ) ( ( Ur ~ U r ) + b ( ~ ) J = 


or 


m, f U,'-U, J + m, f U/-UJ = f n 

( U+'-Uf] = ~f n 


(ring segment) 
(fragment) 


Rotational Impulse-Momentum Law 

mi, [ U'-U, J y l s - U/-UJ a-\> s 
= -% (d-r L )S 


i u ^ 


(ring segment) 
(fragment) 


Kinetic Energy Conservation Law 

"2* Ml, ( U' f + j- ( U* f + J ( Uf J + J- If (ty ) 

= jr m, ( U \) + jr ( (J z /+ j f + f If 


where 


b = ( v - T L ) S 


X = 


»1; 


Ml, * 


(see Fig. 40b) 


(6.19a) 

(6.19b) 

( 6 . 20 ) 


( 6 . 21 ) 

(6.22) 


(6.23) 


(6.24) 


and other quantities retain their previous meanings. 

By using a reduction procedure similar to that described previously and 
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introducing the coefficient of restitution, e, the following equation which 
is analogous to Eq. 6.11 is obtained: 

j [j3 U/-Z rtf ] j = -e[rpU + d UHU f -2<)} (6.25) 

By following a solution procedure similar to that described previously, one 
obtains 


U,'= U, - u 


K l ' ^ 

■+e 


U 2 / = U z - yr\,d ) U 


(6.26) 


(6.27) 


U f — U+ " f w f ' kr iU, 


K u 1 


(6.28) 




2 , I •+ € \ j I 

{ K l > U R 


(6.29) 


where 


HL!2s. + m 6L l ■+■ J3 + "y- 


1 , m, nil * 


(6.30) 


and where all other quantities retain their previous definitions. 


6. 2. 1.3 Governing Equations 

Summarized here, for convenience, are the governing equations of motion 
for both the ring (target) structure and the rotor blade (fragment) . 


Ring-Structure Motion 

As described in Subsection 3.2, the governing equations of motion + for 
either a complete ring or a partial-ring structure may be written as follows 


Rote that Eq. 6.31 represents one of the two forms of the equations of motion 
discussed in Subsection 3.2; the second form, Eq. 3.1, could be discussed in a 
similar fashion if desired. 
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for the spatial finite-element representation of the complete assembled dis- 
cretized structure (CADS) : 

CM] If") + f Pi -[Hjlfl = [ F‘| (6.31) 

where {q*}, {q 1 } represent the generalized displacements and gener- 
alized accelerations, respectively 
[M] is the mass matrix for the CADS 
{P> is an "internal force matrix” which replaces the 
conventional stiffness terms [K] {q} for small 
displacements but also now includes some plastic 
behavior contributions 

[H] represents a "new" stiffness matrix which arises 
because of large deflections and also plastic 
behavior 

{f* } denotes the externally applied generalized forces 
acting on the CADS 

It should be noted that all quantities in Eq. 6.31 refer to the global 
Y,Z inertial reference system indicated in Figs. 41 and 42. For the case in 
which the structure is subjected to distributed linear restoring springs as 
depicted, for example, in Figs. 41 and 42, Eq. 6.31 becomes 

[ M) f f*\ ■+ f Pj + f H] I $1 = { F j “ [ Ks j (6.32) 

where [K s ] represents the global effective stiffness supplied by the elastic 
foundation and/or other "restraining springs". Further, it is presumed that 
Eq. 6.32 has already incorporated within it all pertinent boundary conditions 
and restraints as depicted, for example , in Figs. 42a through 42d. 

As discussed in Subsection 3.3, the timewise solution of Eq. 6.32 may 
be accomplished by employing an appropriate timewise finite-difference scheme 
such as the central difference method. Accordingly, for the cases of CIVM 
fragment impact or of prescribed externally-applied forces, Eq. 6.32 at time 
instant j may be written in the following form: 

CM] ! rij - ( f FT - ffGJ Ifl - 1 Pj - r H ] ! f) ),■ ( 6 . 33 , 
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Let it be assumed that all quantities are known at any given time instant t . . 
Then one may determine the generalized displacement solution at time t (i .e. , 
{q*}j + ^) by the following procedure. First, one employs the timewise central- 
difference expression for the acceleration {q*}_. : 


if). 


I 


( ff*),-.,-*' IV >»*!,,.) 


(6.33a) 


since {q*K is already known from 


1 W 

It follows that one can solve for {q*}_. 

Eq. 6.33 and all other quantities in Eq. 6.33a are known. However, a fragment- 
ring collision may occur between time instants t. and t... ; this would require a 
. . . 3 3 +i 

correction" to the iq*K + ^ found from Eq. 6.33a. Thus, one uses and rewrites 

Eq. 6.33a to form a trial value (overscript T) : 

3” « i • - , .1 

, •+ (*t 

H 


where 


= i "f} ( - + W) i 

Ki,= fft, -i n r , 

= f i\, - 1 f\ 

t i\ } - t A + f* f i 


(6.33b) 


= trial increment 


(6.33c) 


-t- - 




-d "t, = time increment step 

Note that t = j (At) where j * 0, 1, 2, ..., and {Aq*} = 0. Also, no such trial 

J o 

value is needed if only prescribed external forces were applied to the contain- 
ment/deflection ring. 

Let it be assumed that one prescribes at t * t =0 (j«*0) values for the 

# o 

initial velocities {q*} and external forces {F*} , and that the initial stresses 

o o 

and strains are zero. The increment of displacement between time t and time t 

o . 1 

is then given by: 


( 


r i, i .i 


L ( Atl 


(6.34) 


where {q*} can be calculated from 
o 


wherein it is assumed that no ring-fragment collision occurs between t and t 

o 1 

(accordingly, overscript T is not used on {Aq*^ in Eq. 6.34). 
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Fragment Motion 


In the present analysis , the fragment is assumed to be undeformable and, 
for analysis convenience to have semi-circular ends? hence, its equations of 
motion are: 


where 




f 



(Y f ,Z f ) and (Y f , 



0 - 0 


Z f ) denote, respectively, the global coordinates 
and acceleration components of the center of 
gravity of the fragment (see Fig. 43) 

9 represents the angular displacement of the 
fragment. 


(6.35) 

(6.36) 

(6.37) 


In timewise finite-difference form, Eqs. 6.35 through 6.37 become 


v. 

= 

v (6.38) 

T 

= (aZfV 


( A Zf) i 

(6.39) 

T 

= ( * e 



(6.40) 

where overscript "T" signifies 

a trial value which 

requires modification, as 


explained later, if ring- fragment collision occurs between t_. and 

By an inspection procedure to be described shortly, the instant of 
ring-fragment collision is determined, and the resulting collision-induced 
velocities which are imparted to the fragment and to the affected ring segment 
are determined in accordance with the analysis of Subsection 6. 2.1.1 or 6. 2. 1.2. 


6. 2. 1.4 Solution Procedure 

The following procedure indicated in the flow diagram of Fig. 44 (and 
described also in Ref.. 181) , may be employed to predict the motions of the ring 
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and the rigid fragment, their possible collision, the resulting collision- 
imparted velocities experienced by each, and the subsequent motion of each 
body: 


Step 1 


Step 2 


Step ... 3 


Step 4 


Let it be assumed at instant t. that the coordinates {q*}, Y , and 

3 3 

Z , and coordinate increments {Aq*}.,AY , and A Z. are known. One 

3 


j 


j 


can then calculate the strain increments Ac. at all Gauss stations j 

3 

along and through the thickness of the ring from Eq. 4.25. 

Using a suitable constitutive relation for the ring material, the 

stress increments Aa^ at corresponding Gaussian stations can be 

determined from the now-known strain increments Ac. (see Subsections 

3 

2.3 and 4.3.4). Since the a. are known at time instant t. , 

3~l 3“l 

the stresses at t . are given by 0, = c . , + A a . . This information 
3 3 3-1 3 

permits determining all quantities on the right-hand side of 
Eq. 6.33, where for the present CIVM problem {f*K is regarded as 
being zero. 

Solve Eq. 6.33 for the trial ring displacement increments {Aq*}^ + ^. 

Also, use Eqs. 6.38, 6.39, and 6.40 for the trial fragment dis- 

t r r 

placement increments (Ay , (Az^.)^,, and (A0).., 


f'-j+l* v f'j+l' 


j+1* 


Since a ring-fragment collision may have occurred between t^ and 
t^ + ^ , the following sequence of substeps may be employed to de- 
termine whether or not a collision occurred and, if so, to effect 
a correction of the coordinate increments of the affected ring 
segment and of the fragment. 

Step 4a : To check the possibility of a collision between the 

fragment in the vicinity of point A of the fragment 

with ring element i (approximated' as a straight beam) 

as depicted in Fig. 43, compute the trial projection 

(p J of the line from ring node i to point A of the 

, fragment, upon the straight line connecting ring nodes 

i and i-1 , as follows , at time instant t . , , : 

3+1 
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c-M,v, = r Y,- Y A Jj., CM 

f t *1 *t 

l ^ J^'r( 'Jl' 1 )f+l (6.41) 

where the Y,Z are inertial Cartesian coordinates. 

T 

Now, examine <£>^) j +i * three cases are illustrated 

in Fig. 43a. 

T T 

Step 4b ; If (P^) j +1 < 0 or if (p^) > s^^ where s^ > 0, a 
collision between the fragment near point A and 
ring element i is impossible. Proceed to check 
ring element i+1, etc. for the possibility of a 
collision of fragment end A with other ring ele- 
ments . 

T 

Step 4c ; If 0 < (pj j +1 <_ s^, a collision with ring element 

i is possible, and further checking is pursued. 

Next, calculate the fictitious "penetration dis- 
T 

tance" (a.) . of the fragment at end A into ring 
l 3+J. 

element i by (see Fig. 43b) : 

T / T 

( )j*t = l f h,i + ci ( h l£ -h U ) + Jj H - [ d iJ r , (6.42) 


where 


— (hii + a(h 2i - h^^) ] » local semi -thickness of the ring 

element which is approximated as 
a straight beam in this "collision 
calculation" 


“j+l 


C.4c >,v, 




— * tip radius of the fragment at the 
impact end A 


fractional distance of s. from node i 

x 

to where the collision occurs (recall: 


a + 0 


1 , and a. , , should not be con- 
3+1 -r 


fused with the angle (ol ) _. + ^) . 

- tY 16 - 43 ’ 

« the projection of the line connecting 
node i with point A upon a line per- 
pendicular to the line joining nodes 
i and i-1. 
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Step 4d 


Step 4e 


Next/ examine (a^) _. + ^ which is indicated schematically in 
Pig. 43b and is given by Eq. 6.42. 

If (X) ^ $ 0, no collision of the fragment near point A 

upon element i has occurred during the time interval from 

t_. to tj +1 > Hence, one can proceed to check element i+1, 

etc. for the possibility of a collision of fragment end A 

with other ring elements . 
t 

If (a^) > 0, a collision has occurred; corrected coordi- 

nate increments (overscript "C") may be determined approxi-’ 
mately by (see Fig. 43b) : 

C T f “y 

( a Yf ) }+l = ( A V f ) ;+ . -*• (At*) [ Uf- Ufjsin (U-) r , (6 . 44a) 


(AW f )^ r (^t it j[U/-UfJ^S(5 A -)^ l (6.44b) 
0 )j+, = ^ ® + (At J [ u)f - ) (6.44c) 

= fU/-U 2 (6 . 44d) 

^ w,_, ) r , = %•_, ^v.+^f U i~ UJ (6 * 44e) 


(* 5 a - ) r , = ( ^ V K ) ?V( - f U,'~ U, J Sin l d- ) r , (6.44f) 

T 

(* VV X - ) jV| = (-a ** ( At * } [ U- U j C05 ( iv, (6.44g) 


where the after-impact quantities U£, Uj, U^, and may 
be found, respectively, from (Eqs. 6.16, 6.17, 6.13, and 
6.15) or from (Eqs. 6.26, 6.27, 6.28, and 6.29), which- 
ever collision model one wishes to use, and where 


= 


< Ur; ) i 


time interval from actual 

= impact on ring element i 

until t... 

3+1 
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(6.45a) 




points A and C (6.45b) 


The terns, in Eqs. 6.44a through 6.44g, which are 
multiplied by (At*) represent corrections to the 
trial incremental quantities for the (At*) time in- 
terval. Also, since At is small, one may use either 
T 

angle (a.) . . or angle (a.) . in Eqs. 6.44a through 
x 3+1 i 3 

6.44g. 


Step 5 : Having determined the corrected coordinate increments for the im- 
pacted ring element, this time cycle of calculation is now complete . 
One then proceeds to calculate the ring nodal coordinate increments 

and the fragment coordinates for the time step from t. . to t. _, 

3+1 3+2 

starting with Step 1. The process proceeds cyclically thereafter 
for as many time increments as desired. 


If, however, one finds no collision of fragment end A with any 
of the ring segments, the checking process should be repeated 
for any other possible fragment points of impact (such as end B, 
for example) with the ring. 


This solution procedure may be carried out for as many time steps as 
desired or may be terminated by invoking the use of a termination criterion 
such as, for example, the reaching of a critical value of the strain at the 
inner surface or the outer surface of the ring. Appropriate modifications 
of this approximate analysis could be made, if desired, to follow the be- 
havior of the ring and the fragment after the initiation and/or completion of 
local fracturing of the ring has occurred. 


It should be noted that in this approximate calculation, only the coordinate 
increments of the fragment and of the impacted ring segment are corrected. 
Those for all other ring segments are regarded as already being correct. The 
time increment At is regarded as being sufficiently small to make these ap- 
proximations acceptable. 
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6.2.2 Outline of the Collision-Force Method (CFM) 


Although the collision-force method will be described in detail in 
Ref. 183, it is perhaps useful to outline briefly the essence of this approach. 
In this method, the primary information predicted in the collision analysis 
consists of the collision-induced interaction forces themselves. One may 
readily identify the following three versions of the method (where the ac- 
companying sketches indicate the associated qualitative behavior) : 

(a) The deformation of the fragment and the ring may be neglected 
in the collision-interface region. 


Collision 
Starts 
t * t 


o 


During 

Interaction 

V > fc o 




Collision 

Ceases 




where a denotes the total local deformation of the ring and fragment 
at the center of the impact region (for this case a = 0) . 
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(b) The fragment and the ring deform in 
in the collision-interface region. 


a reversible elastic manner 


Collision 

Starts 


During 

Interaction 

t, > t > t 
1 o 


Collision 

Ceases 


t > t. 


A \ 







This reversible elastic local indentation a is related to the 

3/2 

impact force F by the Hertz Law: F <■ ka where k is a con- 
stant found from the elastic properties of the colliding ma- 
terials. 
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(c) The fragment and the ring may deform in a general elastic- 
plastic fashion in the collision-interface region. 


Collision 

Starts 


During 

Interaction 

t > t > t 
1 o 


Collision 

Ceases 



t > t. 




impact 

Force 



Here both the ring and the fragment exhibit permanent deforma- 
tion in the "collision zone"; for this situation the impact 
force may be expressed approximately by F = Na n where N and 
n are constants depending upon the properties of the colliding 
materials . 

For all of the CIVM cases and CFM cases (a) , (b) , and (c ) , the containment 
rings may undergo elastic-plastic membrane and bending behavior. 
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6.3 Application and Evaluation of the CIVM Approach 

Described in this subsection is a set of very limited calculations in 
which the CIVM approach has been employed. The calculations have involved the 
impact of a single rotor blade (fragment) against: (a) a complete free circu- 
lar containment ring and (b) a "fragment-deflection" ring quadrant which is 
supported at a selected station in one of several ways. These brief calcula- 
tions serve to illustrate minimally (1) the effect of the type of collision 
model (consistent mass, CM, or lumped mass, LM) , (2) the effect of the coef- 
ficient of restitution e, (3) the effect of plausible ring material strain- 
rate behavior, (4) calculation convergence with an increasing number of ring 
elements, and (S) preliminary comparisons with experimental response data for 
a complete free ring subjected to impact by a single-blade fragment. These 
studies are described in the following. 

6.3.1 Definition of Example Problems and Calculations 

To illustrate the CIVM approach, the two categories of example problems 
depicted schematically in Fig. 45 were analyzed: (a) a free complete circular ring 
subjected to impact from a single-blade fragment, as represented by the con- 
ditions pertaining to experiment in NAPTC Ring Tests 88 and/or 91 (Ref. 177) 
and (b) a ring quadrant chosen to represent an example deflection device and 
also subjected to single-blade impact. For the complete-ring example, pre- 
liminary experimental ring response and blade motion data (from NAPTC Ring 
Tests 88 and 91, with the latter being the more reliable) are available for com- 
parison with predictions. However, experimental data are not available for the 
ring-quadrant cases. 

Summarized in Table 1 are the pertinent geometric and test-condition 
data for NAPTC Ring Tests 88 and 91. The (complete) free containment ring con- 
sisted of 2024-T4 aluminum; for analysis, its uniaxial static stress-strain be- 
havior was approximated (closely) as being elastic, perfectly-plastic (EL-PP) 
with a yield stress of 50,000 psi and an elastic modulus of 10 psi. A 
single T-58 rotor blade which was fabricated from material designated as SEL-15 
by General Electric was the fragment employed; for the present analysis, this 
"fragment" is treated as being rigid. For the ring-quadrant examples , the same 
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geometric, material property, and test conditions were employed as for the 
Test 88 comple te-ring example. 

Shown in Tables 2 and 3 are the characterizing data for the sequence of 
CIVM calculations carried out for, respectively, (1) the complete ring — 

Runs CR-lB through CR-11B, and (2) the ring Quadrant examples — Runs RQ-1B 
through RQ-9B. For all cases except CR-8B through CR-11B, the conditions 
pertinent to NAPTC Test 88 were used; NAPTC Test 91 conditions were used for 
calculation Runs CR-8B through CR-11B. 

Shown in Table 2 for the complete-ring problem are (a) the number of 
segments into which each quadrant of the ring was discretized for analysis , 

(b) the interpolation function used to represent the displacement field 
throughout each ring element (the linear distribution of the inplane dis- 
placement together with a cubic distribution of the transverse displacement 
is termed LC; the cubic distribution of both the inplane and transverse dis- 
placement along each ring finite element is termed CC for cubic cubic) , (c) 
the type of mass matrix used for each ring element (the mass matrix which is 
variationally consistent with the assumed displacement field is termed C and 
the lumped-mass matrix is termed L) , (d) the ring material property represen- 
tation (elastic, perfectly plastic EL-PP, or EL-PP-SR, where SR denotes 
strain-rate sensitive behavior) , (e) the type of collision model employed 
(consistent mass, CM, or lumped mass, LM) , and (f) the coefficient of impact 
restitution e, where e = 1 represents a perfectly-elastic impact and e = 0 
denotes a perfectly-inelastic impact. 

For the partial-ring or ring-quadrant cases. Table 3 identifies: (a) 
the number of equal- length segments into which the ring quadrant was idealized 
for analysis, (b) the type of mass matrix used for each ring finite element 
(L or C) , (c) the type of support (ideally clamped IC, or smoothly hinged), 
and (d) the type of collision model (CM or 121) used. For all cases, the LC 
type assumed-displacement function over each ring finite element, EL-PP-SR, 
and e = 1 conditions were used in these calculations. 

The various matters examined in these studies are described in the 
following. 
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6.3.2 Effect of Type of Idealized Collision Model 

The two types of idealized collision models, consistent-mass CM and 
lumped-mass LM, discussed in Subsections 6. 2. 1.1 and 6. 2. 1.2, respectively, 
have been employed in example calculations to examine their behavior and rela- 
tive merits. 

First, consider the CM collision model. A comparison of ring quadrant 
calculations RQ-lB vs RQ-2B (see Table 3) have illustrated dramatically the un- 
desirably sensitive nature of the CM collision model. For case RQ-lB, the 
ring quadrant was modeled by 9 equal-length segments while for case RQ-2B, 10 
equal-length segments were employed; in all other respects the conditions were 
the same: the same fragment properties and initial conditions, the ring materi- 
al properties EL-PP-SR, and e = 1 for the coefficient of restitution were em- 
ployed. Dispite this slight change in the modeling of the ring (from 9 to 10 
equal- length elements) , the resulting responses of the ring quadrant and of the 
blade fragment were dramatically different (see Fig. 46) when the CM model was 
employed. Increasing the number of segments to 15 to model the ring quadrant 
in case RQ-3B resulted in ring quadrant and blade responses that also differed 
significantly from the 9-element case but very little from the 10-element case. 
These results illustrate the sensitive nature of the CM collision model. 

On the other hand, the use of the LM collision model for identical ring- 
quadrant modelings: 9, 10, and 15 elements for cases RQ-4B, RQ-5B, and RQ-6B , 
respectively, demonstrated improved convergence and much less sensitivity to 
the number of segments (elements) used to model the ring quadrant. These re- 
sults are illustrated in Fig. 47. Note that the 9-element result differs some- 
what from the 10-element result, but the 10-element result and the 15-element 
result are nearly the same; this suggests that the 10-element calculation pro- 
vides a "converged" solution. Because of the greater consistency of the LM 
collision model in comparison with the CM collision model, the LM collision 
model has been selected for principal use, although additional CM collision 
calculations have also been carried out. 

Similar predicted ring response and blade response behavior for single- 
blade impact upon a free complete ring is illustrated in Fig. 48. CM collision 
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model cases CR-1B and CR-3B for ring discretizations 5-5-5-5 and 9-6-6-6 , 
respectively, are compared with LM collision model case CR-5B which utilizes 
ring discretization 10-6-6-6 (see Table 2) ; the latter calculation is assumed 
to provide a converged result. Here again the CM results appear to exhibit 
somewhat "erratic" behavior. 

An approximate collision model which leads to very different results 
when only a slight change in modeling of the ring is employed is undesirable 
and unacceptable . A trustworthy approximate collision model should not exhibit 
such extreme sensitivity to a slight change in the discretization modeling of 
a containment ring or of a "deflector” ring structure. Accordingly, the CM 
collision model has been set aside and the LM collision model has been adopted 
for future calculations until a better model has been devised. 


The reasons for the erratic CM collision model behavior versus the more 
satisfactory LM collision model features may be seen most conveniently perhaps 
by examining the collision-imparted velocity behavior of the blade and the 
affected ring element. It should ne recalled that the ring element was 
idealized by a straight-line beam element (for the impact-collision calcula- 
tion only) . Consider the case in which the uniform beam element with length s , 
mass m^, is struck by the blade tip at a point very near one end of the beam 
element (for example, near node 1, see Fig. 40) , in which case a = 0 and 8 = 1; 

assuming that the normal impulse resulting from the collision is p , the blade 

n 

c.g. normal-direction velocity change and angular velocity change due to impact 
are given by Eqs. 6.2 and 6.4 (also Eqs. 6.20 and 6.22) as 





Ir 


(6.46a) 


(6.46b) 


where the angle 6 is the inclination of the blade to the normal of the beam 
element. 


The beam nodal velocity awl the angular velocity changes predicted by 
the CM collision model (Eqs. 6.14, 6.16, and 6.17) for a collision occurring at 
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node i(a = 0 and 3 = 1) are 


y an 1 1 A. _ 4 1 




m. 


. , i 7 n ] _ f , , y™ _ g i ft, 

( w / cM I J + yam 2 j KKl r w r 


_ ± _k 

CM J, +yam MV 5 n V 


( u) t '-(A)J= —2- 


(6.47a) 


(6.47b) 


(6.47c) 


On the other hand, the LM collision model predicts the beam nodal velocity and 
angular velocity changes (Eqs. 6.26, 6.27) to be 


( (J, U, J LjV1 
( U/- U 2 ) lM 

( cO r — cOjr ) L/l/ f 


= P-k, r = 

= = 
m r 

I ft, 



0 


(6.48a) 


(6.48b) 


(6.48c) 


Comparing Eq. 6.47c with Eq. 6.48c, it is seen that, for a certain normal im- 
pulse p^, the angular velocity change predicted by the CM model would be six 
times larger than that predicted by the LM model. This larger consequent 
change of the orientation of the beam segment for the CM model compared with 
that for the LM model may subsequently lead to the larger change of the angle 
0, and hence, the angular velocity change of the blade (see Eq. 6.46b) for the 
subsequent impact (s) of the blade with the beam element will also be more 
seriously altered. As observed in Fig. 46, the use of the CM collision model 
with 10-element and 15-element modeling of the ring quadrant, after a certain 
time stage, predicts that the angular velocity of the blade becomes negative 


157 



(the direction of rotation is reversed) . No such reversal of the direction of 
blade rotation is predicted when the LM collision model is used in corresponding 
10-element and 15-element calculations for the quadrant ring. 

The above example case may serve to illustrate (a) the relative insensi- 
tivity of the LM vs the CM idealized collision analysis and (b) associated 
sources of the differences in predicted collision behavior. 

Finally, it should be noted that at initial impact: (a) the angle 6 be- 
tween the blade and the normal to the impacted ring segment and (b) the loca- 
tion of the impact point on that impacted straight-line collision-model ring 
segment (as reflected by the associated values of a and 3) will be different 
when one uses different numbers of ring segments to represent the ring. As 
the number of ring segments increases, 0 will approach the true value for the 
ring itself, but a and 3 will continue to "fluctuate”. These factors to- 
gether with the sensitivity of the CM collision model are responsible for the 
significant difference in predicted blade motion between the 9-element CM 
case RQ-1B and the 10-element CM case RQ-2B. 

6.3.3 Effect of Coefficient of Restitution 

As noted in Subsection 6. 2. 1.1, the coefficient of restitution e is a 
useful quantity for discussing "local impact” of two colliding bodies, since 
this quantity represents the ratio of the relative velocity of the two bodies 
after impact to that before impact. Perfectly elastic impact is represented 
by e = 1 whereas e = 0 denotes perfectly-inelastic impa'ct (i.e., zero relative 
velocity after impact) . However , the importance of the value of e upon the 
transient structural deformations of a fragment-impacted structure such as the 
containment ring depicted in Fig. 45a is not obvious. 

To examine this matter, CIVM calculation examples have been carried out 
by using e «* 1 and e = 0 in (1) cases CR-1B and CR-2B, respectively, (2) cases 
CR-5B and CR-6B, respectively, and (3) cases CR-10B and CR-llB, respectively. 
Although these three categories of comparisons involved somewhat different 
modelings, as can be seen from Table 2, it was found that to the scale of the 
plots shown on Fig. 48, the differences in ring structural responses for e = 1 
vs e = 0 in each category are almost imperceptible . A more critical comparison. 
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however, is afforded by examining the predicted strains at given locations; 
shown in Fig. 49 are strains predicted on both the outer surface and the inner 
surface of the ring at two interesting 0-locations (0 = 13.5° and 0 = 49.5°) 
for case CR-10B (e * 1) and for case CR-11B (e = 0) — these cases involve the 
best combination of modeling conditions of all cases listed in Table 2. It is 
seen that generally the circumferential strains predicted for e - 1 are slightly 
larger than when e =* 0 is employed. It should be noted that initial impact 
occurs at 0 = 44.5°; in regions or quadrants well removed from the impact zone 
or quadrant, the predicted strains are very small. For example, at 6 * 233°, 
the peak outer surface strain is only about 3.35 x 10 3 in/in. 

Another interesting effect of the coefficient of restitution is that the 
number of ring-fragment collisions which has occurred up to a given instant 
in time after initial impact is much larger for perfectly inelastic impact 
(e = 0) than for perfectly elastic impact (e = 1) . This is illustrated in 
Fig. 50 where the accumulated number of impacts (ANI) is shown as a function 
of time after initial impact for cases CR-10B and CR-11B. For perfectly-elastic 
impact, e = 1 (case CR-10B) , the ANI reached 73 by 370 microseconds after initial 
impact; no further ring-fragment impacts occurred thereafter (at least during 
the period to 818 psec after initial impact, during which the response was ex- 
amined in these calculations). On the other hand, for perfectly-inelastic im- 
pact, e = 0 (case CR-11B) , the ANI reached 377 by 392 microseconds after initial 
impact, with no further impacts occurring to at least 818 psec after initial 
impact. 

The large ANI for the perfectly-inelastic case (e = 0) as compared with 
the perfectly-elastic collision calculation (e = 1) may be readily appreciated 
by recalling that for e ■> 0, the relative velocity of the contact points of the 
blade and the impacted ring segment are zero. For the next time increment. At, 
of the prediction process, the ring and the blade are treated in a trial incre- 
mental calculation as moving independently of each other. Thus, the motion of 
the impacted ring segment is "retarded" by the action of the internal forces 
applied to it by its neighboring ring elements, while the fragment proceeds 
with "uninhibited" motion. The subsequent collision inspection, therefore, 
frequently reports that during this At another collision has taken place — and 
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the process continues. But for e = 1, a nonzero relative velocity of the im- 
pact points of these two objects occurs; subsequent collisions, hence, are less 
frequent. 

It is interesting also to note in the following tabulation that after the 


impacts have ceased, the blade has somewhat different cg-translational and ro- 
tational velocities for the e = 1 case as compared with the e - 0 calculation: 


Blade Motion Quantities 

Pre- Impact 
Conditions 
of Blade 

After Impacts have Ceased 

e = X 

Case CR-10B 

e = 0 

Case CR-11B 

(in/sec) 

- 7884 

- 2826.8 

- 2138.4 

V (in/sec) 
z 

0 

- 3038.2 

- 3015.4 

Q)^( rad/sec) 

+ 1638.3 

+ 3137.5 

+ 3015.4 

a f (deg) 

0 

+ 111.4* 

+ 105.9** 

a^(deg) at TAII = 818 Msec 

- 

+ 192.0 

+ 170.3 


It is seen that after all impacts have ceased, the blade has larger cg- 

translational (V and V ) and rotational (us.) velocities for the e = 1 case 
y 2 i 

than for the e = 0 case. Finally, had these calculations been carried out for 
longer times, it is clear from the fragment and ring motions that further 
blade-ring collisions would have been seen. 

6.3.4 Strain Rate Effect 

In the present analysis the fragment (blade) is treated as being rigid, 
but the containment ring is deformable. Thus strain-rate dependent mechanical 
behavior of the ring can be taken into account, and is expected to influence 
primarily the response of the ring and secondarily the motion of the fragment. 

As described in Subsection 2.3.3, it is assumed herein that strain rate 
E ■ raises the uniaxial yield point of the material, approximately as follows: 

— 

At time after initial impact (TAII) = 370 Msec, or time = 552 + 370 = 922 Msec 
At TAII = 392 psec 
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where 


O 

o 




strain-rate dependent uniaxial yield stress 
static uniaxial yield stress 


and D and p are constants which depend upon the material involved. For 6061-T6 
aluminum, D = 6500 sec ^ and p = 4 are commonly used; these values have been 
employed in the cases denoted by SR (in EL-PP-SR) in Tables 2 and 3. It is also 
assumed that the elastic modulus of the material at room temperature conditions 
is not affected by e. 

The effect of £ in the ring-blade problem can be seen by comparing the 
predicted ring and blade responses for case CR-5B (EL-PP) versus case CR-7B 
(EL-PP-SR) ; these cases are identical except for the inclusion of the strain- 
rate effect in case CR-7B. Let it suffice for present purposes to note that 
the inclusion of the £ effect in case CR-7B manifests itself in making the 
ring "stiffer" than in case CR-5B (EL-PP) ; the peak deformations of the ring 
were reduced slightly by including the £ effect. The result is in accord with 
the more detailed e-affected structural responses described in detail in Sub- 
section 5. 3.6. 3. 


Finally, some blade-motion data (which is of secondary interest) for 
cases CR-5B and CR-7B are compared in Table 4; these include the fragment cg- 
velocity components and and the angular velocity U)^ of the fragment 
(blade) . From Table 4 it is seen that the predicted velocity histories of the 
blade for the EL-PP and the EL-PP-SR calculations are very similar. In fact, 
the "residual" velocities of the blade after no further impacts occur (up to 
at least TAII * 818 Msec) differ only slightly are are as follows: 
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Case 

TAII for 
Last Impact 
(Usee) 

AMI 

(in/sec) 

V (in/sec) 
z 

(rad/sec) 

CR-5B 
EL-PP 
e = 1 

348 

89 

- 2791.0 

- 2999.4 

+ 2611.7 

CR-7B 
EL-PP-SR 
e = 1 

327 

47 

- 2912.0 

- 2991.3 

+ 2621.5 


The most interesting difference in these results is that the accumulated number 
of impacts to "impact cutoff" is 89 and 47 for the EL-PP case and the EL-PP-SR 
case, respectively. That is, the inclusion of "strain-rate stiffening" of the 
ring reduced the predicted number of impacts (up to at least TAXI = 818 usee) ; 
this, incidentally, reduces the computational effort and time. 

6.3.5 Comparison of Predictions with Experiment for a 

Complete Free Ring Subjected to Rotor Blade Impact 

Having examined the effects of various features in the CIVM approach 
and various modelings of a free complete circular "containment ring" subjected 
to impact by a "nondeformable" single rotor blade, it is now perhaps of in- 
terest to compare predictions with measurements from NAPTC Tests 88 and 91 
(see Table 1) . These comparisons are regarded as tentative because of the 
facts that: 

(a) the time from the instant of initial impact until the 
"time instant" of each photograph in each experiment 
is uncertain by as much as, perhaps, about 200 to 300 
microseconds* (that is, initial impact may have occurred 
earlier than reported in Ref. 184 by about 200 to 300 
Usee) and 

(b) the present CIVM analysis contains some simplifying 
assumptions which expedite the present study but 
which do not accommodate some aspects of the antici- 
pated and observed behavior (these deficiencies could 
be remedied in future work) . 

_ 

This uncertainty will be essentially eliminated in planned forthcoming tests 
of this type through the use of improved techniques. 
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The predictions discussed in this subsection are limited to only the 
best two calculation cases for each experiment. Tests 88 and 91. Accordingly, 
shown in Fig. 51 are deformed-ring profiles and blade locations at a sequence 
of times for : experiment (NAPTC Test 88} , case CR-5B (EL-PP , e = 1) , and case 
CR-7B (EL-PP-SR, e = 1) . The plotted points represent the locations of the 
node points or ends of the finite elements or segments into which the ring 
has been discretized for analysis. For comparing predictions with experiment, 
the initial impact location was matched to that observed experimentally. It 
is seen that there is reasonably good agreement between predictions and ex- 
periment, with calculation case CR-7B providing the better comparison. In 
Fig. 51, the time interval between the "instant of impact" and the photograph 
for which the deformed-ring data were measured is as reported in Ref. 184. 

Figure 52 shows a similar comparison between predictions and experimental 
results measured for NAPTC Test 91s lie predictions are for case CR-10B (EL-PP-SR, 
e = 1) and case CR-llB (EL-PP-SR, e = 0) . Hence, seen again are the effects of 
perfectly-elastic impacts (e = 1) versus perfectly-inelastic impacts (e = 0) ; 
on the scale of Fig. 52, this coefficient-of -restitution effect is minor. Here 
also (1) it is seen that the predictions agree favorably with experiment , and 
(2) the "experimental instant of impact" used is as reported in Ref. 184. 

A preliminary assessment in Ref. 185 of the Ref. 184 data for NAPTC 
Tests 88 and 91 suggests that the actual instant of impact in each case may 
have been significantly earlier than reported in Ref. 184. For NAPTC Test 91, 
for example, it was estimated in Ref. 185 that the actual instant of impact was 
perhaps about 200 to 300 microseconds earlier than implied in the experimental 
data of Fig. 52 and Ref. 184. Hence, as a matter of curiosity, an "impact time 
correction" of 240 psec has been applied to the NAPTC Test 91 data and the re- 
sulting comparisons with the case CR-1QB predictions are shown in Fig. 53. It 
is seen that the predicted ring responses are somewhat larger than the Test 91 
data show when a revised' instant of initial impact is used. As mentioned in 
the beginning of this subsection, this discrepancy may possibly be due to the 
uncertainty of the time elapsed from the instant of actual initial impact to 
the "test-recorded" initial impact instant, and may also be due to the simpli- 
fying assumptions employed in the present approximate predictions. Further 
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experimental and theoretical studies of this situation are required. 

With reference to the results shown in Figs. 51-53, it is seen that the 
predicted ring responses agree fairly, well with experiment, but the observed 
blade motion (which is of secondary interest in the present context, but of 
considerable interest in other situations) is not in satisfactory agreement 
with experimental observations. As for the ring itself, the method of large 
deflection, elastic-plastic transient response analysis has been evaluated 
and its accuracy verified by comparison with reliable experimental data in 
Subsection 5.3 and, hence, is not a source of significant error. However, 
one source of error is readily apparent: in the experiments, the rotor blade 
undergoes a significant amount of deformation over a portion of its length at 
and near the "impact tip" during a brief period following initial impact; 
little further blade deformation is observed at later times. The blade, hence, 
becomes a shortened fragment with a smooth curved portion at one end while the 
other end essentially retains its pre-impact configuration. Accordingly, the 
"new" fragment which continues the blade-ring interaction process often has 
the same mass as before but has a reduced moment of inertia — this factor 
probably accounts for much of the discrepancy between the observed and the pre- 
dicted blade motion. Other neglected factors which may contribute to the ex- 
perimental-theoretical discrepancies include: (a) the neglect of the true ma- 
terial properties of the blade in the collision calculation itself (but this 
should be minor as results from the e = 1 and e = 0 extremes show) and (b) 
forces arising with and energy loss by gouging the ring by the blade. 

6.3.6 Responses of a Variously-Supported Ring 
Quadrant to Rotor Blade Impact 

As noted in Subsection 6.3.1 and in Figs. 36b and 45b, partial rings sure 
of interest as possible fragment deflection devices. Accordingly, illustrative 
calculations have been carried out for certain ring quadrant configurations 
(see Fig. 45b) subjected to impact by a single rotor blade. These ring-quadrant 
configurations are supported in one of the three following ways: (1) ideally 
clamped (IC) at 0 = 0° , (2) smoothly hinged (H) at & = 0°, or (3) smoothly 
hinged (H) at 9 = 27°. The pertinent CXVM calculation features are suraaarized 
in Table 3, and include cases RQ-1B through RQ-9B; identical "blade release" 
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conditions at t = 0 are used for all cases. The Test 88 conditions reported in 
Table 1 were used for these calculations. No experimental data, however, are 
available for comparison with these predictions. 

For the ideally-clamped quadrant ring, the predicted ring quadrant and 
blade responses have already been noted for cases HQ- IB through RQ-3B and cases 
KQ-4B through FQ-6B in Figs. 46 and 47, respectively. The results from case 
RQ-7B differ little from those of case RQ-6B and hence are not shown here. 

Compared in Fig. 54 are the responses of quadrant rings with an ideally- 
clamped end (case KQ-5B) and a smoothly-hinged end (case RQ-8B) . These responses 
exhibit the expected qualitative differences. The "blade diversion angle" 

(<j>, see Fig. 36b) is predicted to be approximately 42.6 s pox both the IC case 
RQ-5B and the hinged-end case RQ-8B. 

The influence of locating the smoothly-supported hinge at 0 * 27° 

(case RQ-9B) rather than at the end, 0 = 0° , of the quadrant ring is depicted 
in Fig. 55 by the sequence of deformed ring quadrant configurations for case 
RQ-9B. For this case, the rotor blade cg-trajectory diversion angle <J> is 
43.3°. 

These examples, RQ-1B through RQ-9B, illustrate only a few of the many pos- 
sible situations (various ring boundary conditions, various fragments, etc.) and 
configurations (variable thickness ring, etc.) which may be worthy of study for frag- 
ment deflection purposes, and to which the CIVM analysis can readily be applied. 

6.4 Comments on Potential Extensions of the CIVM 
Approach to More Complex Problems 

In this report the CIVM approach has been illustrated and applied to only 
a few (see Fig. 45) of the many situations of interest; also, some simplifying 
approximations have been employed. One fragment having a simple shape has 
been used, rather than (a) many fragments (Fig. 37c) or (b) a fragment of 
complex geometry (Fig. 37b) ; the useful but drastic simplifying assumption that 
the fragment acts as a nondeformable body was employed . At this point it is 
useful perhaps to comment briefly on steps which might be taken within the CIVM 
context to perform a more physically-realistic analysis and/or to analyze con- 
tainment/deflection devices which are subjected to more complicated fragments 
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and/or conditions. 

It is convenient for discussion to consider two categories of problems: 

(1) those in which various kinds and numbers of fragments are involved but in 
which the fragments are all nondef ormabl e and (2) a category analogous to (1) 
except that all of the fragments are treated as deformable . There could also 
be analyses in which for a certain period of time in the impact-interaction- 
response behavior, the fragment (s) is regarded as deformable, and in another 
appropriately-selected time period as nondef ormable . In addition to consider- 
ing problems in which all motion and deformations occur in one plane (planar 
problems) as is the case throughout this report, more general motion and de- 
formation conditions could be analyzed. However, discussion herein shall be 
restricted to planar problems. 

First, let attention be directed to cases involving only nondef ormable 
fragments. Beyond the simple single- fragment problems depicted in Fig. 36, 
one could apply the CIVM approach readily to analyze the interaction with con- 
tainment-deflection rings of complex single fragments such as those shown in 
Fig. 37b or multiple fragments such as shown in Fig. 37c. Only the bookkeeping 
becomes more complicated and lengthy than for the Fig. 36 category of problems; 
no conceptual or analytical difficulties are anticipated. If in addition to 
interacting with the ring (and/or with each other) the fragment(s) collide 
with one or more blades of the "remaining" rotating rotor and hence receive a 
new kick, the bookkeeping becomes still more lengthy and coup lex, but no basic 
difficulty is anticipated . However, as noted earlier, treating the fragment (s) 
as nondeformable will lead to progressively less accurate blade-fragment motion 
and accordingly to less reliable predictions of subsequent blade-rotor collisions, 
blade-ring collisions, and ring structural response. How important this effect 
will be upon the predicted containment/deflection capability of a given system 
cannot be ascertained at present. However, it is believed that some modification 
of the analysis should be made to provide more accurate blade-motion predictions 
for the time period following the initial collision-interaction of the fragment 
with the containment/deflection ring than is currently afforded by the present 
rigid- fragment analysis. 

Next, therefore, it is appropriate to consider modifying the analysis 
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to account for deformable fragments. Starting with the simple single-fragment 
configurations of Fig. 36, one could proceed logically to treat the blade as 
deformable by modeling it as consisting of an appropriate number of finite 
elements (or segments) . While somewhat more complex inspection procedures 
would need to be devised to determine when and where ring-fragment collisions 
would occur, the same basic collision-interaction calculation as now employed 
could be used; as done for the ring itself at present, the blade could be 
treated as behaving in an EL-SH-SR manner. In this scheme finite-element 
modeling of both the ring and the blade could be employed throughout the analy- 
sis from initial impact until all impact/interaction behavior of interest had 
occurred; more degrees of freedom and more computing time them currently needed 
in the present rigid-blade analysis would be involved. Some economies could 
perhaps be effected by treating the blade as deformable in the initial impact 
stage and at other critical collision-interaction stages, and by treating the 
blade as nondeformable for the "in-between" periods. 

Another useful variant of this analysis extension might be to treat the 
ring always and the fragments (always or "intermittently") as deformable but 
to treat the still-spinning rotor disk and the remaining attached rotor blades 
as nondeformable. Further study is needed to determine an optimum analysis mix 
for both accuracy and computational efficiency. 

Finally, cognizance should be taken of one aspect of the immediate col- 
lision analysis which should be improved — the seriousness of this defect with 
respect to predicting efficient, assured containment is not entirely clear. 

In the present study, the fragment-ring collision has been idealized as a process 
in which the applied and reactive forces could occur only in a direction per- 
pendicular to the surface of the impacted ring segment. For certain cases, this 
may be a very good approximation; for other cases (where gouging of the ring by 
the blade is imminent or present, for example) , there would be also both applied 
and reactive force components in a direction parallel to the surface of the ring. 
Hence, a modification of the collision analysis to predict and to take into 
account such "shearing forces" — as in machining and tool-wear studies — would 
be a useful step toward a more realistic simulation of the actual physical situ- 
ation for certain cases. 
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SECTION 7 


SUMMARY AND CONCLUSIONS 


7 . I Summary 

The present study is devoted principally to developing and validating a 
spatial finite-element (FE) variational method for analyzing the large-deflection 
elastic-plastic transient deformations of simple structures. As a result, ac- 
curate FE predictions of transient strains and large transient deformations of 
simple structures subjected to known forcing functions have been demonstrated 
(see Section 5 and Subsection 7.2). A practical problem to which the present 
method of analysis has been applied is that of containment/deflection struc- 
tural ring responses to engine rotor-blade fragment impact. 

The equations which govern the large deflection dynamic and/or static 
responses of a solid continuum are discussed in tensor form for convenience and 
generality. Elastic, elastic-plastic, and strain-rate dependent material be- 
havior are considered. Under the restriction that the strain is small, con- 
siderable simplification is attained by assuming that the strain can be de- 
composed into an elastic (thermodynamically reversible) part and a plastic 
(thermodynamically irreversible) part. The Mises-Hencky yield criterion and 
its associated flow rule are adopted to describe the elastic-plastic behavior 
of an initially homogeneous isotropic material. The strain-hardening behavior 
is taken into account by using the mechanical sublayer model. The strain-rate 
effect is approximated by assuming that the uniaxial stress-strain relation is 
affected by strain rate, only by a quasi-steady increase in the yield stress 
above the static-test yield stress. 

The spatial assumed-displacement finite-element (FE) approach is used to 
approximate the true infinite degree-of-freedom description of the actual con- 
tinuum by one which involves a finite number of degrees of freedom. The finite- 
element concept is used in conjunction with the Principle of Virtual Work and 
D'Alembert's Principle to obtain the equations of motion of a general solid con- 
tinuum which is permitted to undergo large-deflection elastic-plastic transient 
deformations. The resulting equations of motion are developed in two forms: 

(a) the conventional form, and (b) an improved form; the latter represents a 
new development. In both forms, the Lagrangian description for displacements, 
strains, and stresses is employed; that is, the initial undeformed configuration 
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of the continuum is used as a reference for the subsequent motion of the continu- 
um. The resulting equations of motion consist of a finite-size system of second 
order ordinary (coupled) nonlinear differential equations with the unknowns to 
be determined being the values of the degrees of freedom (generalized displace- 
ments) at the nodes of the finite-element assemblage which represents the con- 
tinuum. This set of equations is solved timewise by using a direct numerical 
integration scheme with an appropriate temporal finite-difference approximation 
operator. 

The formulation is developed in detail for a curved beamlike structure 
which undergoes planar deformations with: (a) zero transverse shear deformation 
(Bernoulli-Euler-type of deformation behavior) or (b) nonzero transverse shear 
deformation (Timoshenko- type of deformation behavior) . The nonlinearities re- 
sulting from both large-deflections and elastic-plastic material behavior are 
included. 

For the Bernoulli-Euler-type of curved beam element, two sets of assumed 
displacement functions have been employed: (1) both the axial and the transverse 
displacements are represented by cubic order polynomials with small-amplitude 
rigid-body modes included (this is termed a CC element) and (2) the axial dis- 
placement is represented by a first-order polynomial while the transverse dis- 
placement is represented by a cubic-order polynomial, with small-amplitude 
rigid-body modes included (this is termed an LC element) . 

As for the Timoshenko- type of curved beam elements , four sets of assumed 
displacement functions each with various orders of polynomials to approximate 
the transverse displacement have been studied. 

An assessment of this method of analysis is made by means of a sequence 
of problems for beam and ring example structures which are subjected to transient 
mechanical loading or to initial impulsive loading; the present predictions are 
compared with reliable experimental data and/or independent predictions (finite- 
difference and/or analytical) . The temporal central-difference finite-differ- 
ence operator is used for most of the calculations; however, the use of Houbolt's 
method and Newmark 1 s method is also explored for the present nonlinear problems. 

The structural response predictions and evaluations were carried out in 
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the following categories: 


(a) Small-Deflection Linear-Elastic Transient Response 

(1) An impulsively-loaded simply-supported beam was analyzed 
by the present FE method with the Bernoulli -Euler type of 
LC element and the predicted responses were in excellent 
agreement with the exact normal mode solution. This com- 
parison served principally to verify the correctness of 
the computer program. 

(2) A free-free "thick" beam loaded by a transient concentrated 
load at midspan was analyzed by using the several types of 
Timoshenko-type finite elements developed in this study. 

The present predictions were compared with the available 
exact modal solution. 

(b) Large-Deflection Elastic-Plastic Transient Response 

Included in this category were the following impulsively- loaded 
structures: (1) a clamped- clamped beam, (2) a ring sector with 
clamped ends, and (3) a free ring. Transient response pre- 
dictions carried out by using the present finite-element ap- 
proach were compared with experimental measurements (of 
transient strains and large transient def ormations ) and 
finite difference predictions . 

The problem of engine rotor fragments interacting with either a complete 
(containment) or a partial (deflection) structural ring is discussed as an ex- 
ample application of the present analysis to a problem of current practical in- 
terest. Energy and momentum considerations are employed to predict the collision- 
induced velocities which are imparted to the colliding fragment and to the af- 
fected ring segment; the associated analysis method is termed the collision- 
imparted velocity method, CIVM. This collision analysis is combined with the FE 
analysis developed in this study to permit one to predict the resulting large 
deformation responses of containment/def lection rings. Comparisons with limited 
experimental data are also given. 
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7 .2 Conclusions 


A spatial assumed-displacement finite-element variational formulation 
and solution scheme for the prediction of large-deflection elastic-plastic 
transient responses of curved beamlike structures has been developed and has 
been implemented in a computer program; both the older conventional FE formu- 
lation and a new improved FE formulation for analyzing this class of problems 
have been developed and discussed. The accuracy and versatility of these FE 
methods of analysis have been demonstrated on problems for beam and ring ex- 
ample structures which are subjected to transient mechanical loading or to 
initial impulsive loading. The present finite-element-predicted responses are 
found to be in reasonably good agreement with experimental measurements and/or 
with independent predictions (finite-difference and/or analytical) . Also, the 
present finite-element analysis combined with the impact analysis has been 
applied to the prediction of containment/deflection structural ring responses 
to engine rotor fragment impact; the predicted ring responses agree favorably 
with experiment. 

On the basis of the present study, the following conclusions may be 

stated: 

(1) The improved finite-element formulation is more efficient and 
economical for computing the large-deflection elastic-plastic 
transient responses of simple structures than the conventional 
finite-element formulation. 

(2) The 3-point central-difference time integration method gives 
very accurate amplitude and phase transient predictions as 
long as the time increment size used is small enough. The 
largest permissible time increment size At which will avoid 
calculation instability is affected by the severity of the 
structural responses; that is, large-deflections play the 
key role in stiffening the structure and thus requiring a 
smaller At than that required for small deflection transient 
response problems to avoid numerical instability. Also, the 
large-deflection effects render Houbolt's method and Newmark's 
method no longer "unconditionally stable" as they are for 
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small-deflection problems. This is consistent with the 
observations in Ref. 75. 

(3) The Bernoulli -Euler-type of curved beam element with CC 
assumed-displacement functions exhibits significantly im- 
proved predictions, especially for the strain, over the 
Bernoulli-Euler-type of curved beam element with I C 
assumed-displacement functions for a given number of de- 
grees of freedom. Also, the former converges very rapidly, 
but at the expense of a smaller allowable time increment 
step size compared with the latter. 

(4) The Timoshenko-type of beam element with linearly-varying 
assumed-displacement functions (Tl) can provide accurate 
small-displacement linear-elastic transient response pre- 
dictions only if the element size is kept small enough. 

However, in order to obtain more accurate coarse-mesh 
solutions, one would need to employ assumed displacement poly- 
nomial functions of higher order. Better and more efficient solu- 
tions, however, are obtained by using T2, T3, or T4 elements, with 
the T4 type of element being the best. Unfortunately, however, 
pertinent experimental and/or predicted results for large- 
deflection, elastic-plastic transient responses with im- 
portant transverse shear deformation behavior against which 

to compare predictions which could be readily accomplished 
with the present FE analysis have not been located. 

(5) It is substantially more efficient to use the lumped-mass 
matrix version (L) rather than the consistent-mass matrix 
version (C) of the finite-element method for typical tran- 
sient response problems of the present class; the allowable 
At is larger for the former than for the latter. Further, 
since the transient responses predicted by using the C and 
the L finite-element calculation are quite close to each 
other, it is recommended that the lumped-mass (L) version 
of the finite-element method be selected as being more 
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efficient (and adequate) for engineering prediction purposes. 

(6) For the present planar structures, it has been demonstrated 
that the present FE method provides accurate predictions 

of transient structural responses , both deformations and 
strains, for strains at least as large as about 5 per cent. 

(It is reasonable to hope that reasonably accurate strain pre- 
dictions will be provided by this method for strains as large 
perhaps as 20 per cent which approaches fracture levels of 
strain for many common aerospace structural materials.) 

(7) Based upon the limited comparative studies reported in 
Subsection 5. 3. 7. 3 between the finite element and the 
finite difference method, the lumped-mass improved- 
formulation finite-element method is competitive with 
the finite difference method with respect to efficiency 
(confuting time and cost) and accuracy for predicting 
large-deflection elastic-plastic transient responses of 
simple structures. 

(8) The use of the present approximate CIVM approach 
(Section 6) , in conjunction with either the present finite- 
element procedure or the finite-difference method provides 
reasonable estimates of the fragment and ring responses 
arising from single-blade-fragment impact upon a struc- 
tural ring. Several deficiencies in the present "first- 
cut" CIVM analysis could be remedied readily by further 
work, and would clearly lead to improved qualitative and 
quantitative predictions. 

7.3 Suggestions for Further Study 

Possible areas of interesting and useful further research along the 
lines of the present investigation cure described in the following. 

Although one can always resort to numerical experimentation to provide 
a suitable time increment step size to insure stability, it would be desirable 
to develop criteria for the maximum allowable time increment step for the large- 
deflection nonlinear structural transient response problems of the present type; 
a useful beginning for this type of analysis is described in Chapter 3 of 
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Ref. 186. Also, before more definite conclusions can be drawn as to whether any 
one timewise operator is superior to the others, more extensive studies than 
could be carried out here are required. 

Pertinent theoretical and/or experimental results for large-deflection 
elastic-plastic transient responses of beamlike structure having significantly 
large transverse-shear deformations have not been found for checking the present 
Timoshenko- type of finite element; the evaluation of the present Timoshenko- 
type element for this kind of problem, therefore, remains to be accomplished. 

The present-developed finite-element analysis method for predicting 
large-deflection elastic-plastic transient structural responses is quite gen- 
eral and pertains to any type of loaded body. However, in this report, its 
application is demonstrated only for simple (curved and/or straight) beamlike 
structures whose significant deformations are confined to one plane. An ex- 
tension and application of the present method to analyze plate and shell struc- 
tures or to structures with complicated geometric shape , material properties , 
and boundary conditions would be of considerable interest and value, and would 
provide useful versatility to permit analyzing this group of nonlinear three- 
dimensional-deformation transient-response problems — which comprises the 
largest group of practical-interest problems of this type. Although there 
are many ways of constructing a lumped-mass matrix (fro-, finite-difference 
equations, for example) there needs to be developed for these general dynamic 
systems, better rational ways of constructing the lumped mass matrix. 

More definite conclusions need to be drawn as to which of the two gen- 
eral methods the finite-difference or the finite-element method will prove to 
be superior for particular types of problems in this nonlinear transient re- 
sponse category. 

Further, the combining of the finite-difference and the finite-element 
procedures in space to take maximum advantage of the special merits of each 
method for appropriate parts of the structure (that is, one might use the 
finite-difference procedure for smoothly-varying regions of the structure and 
the finite-element procedure in regions of structural irregularities) is also 
suggested for further research. 
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Finally, as noted in Subsection 6.4, the extension of the present 
collision-imparted velocity method to perform a more physically-realistic 
analysis (for example, to account for the roughness of the ring-fragment 
inpact surfaces, the deformations of the fragment, etc., and to analyze con- 
tainment/deflection devices which are subjected to more complicated fragment 
impact and/or conditions) would be a useful variant of this collision- imparted- 
velocity method of analysis. 
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TABLE 1 

DATA CHARACTERIZING NAPTC RING TESTS 88 and 91 



Test 88 

Test 91 

Outside Diameter (in) 

17.619 

17.619 

Radial Thickness (in) 

0.153 

0.152 

Axial Length (in) 

1.506 

1.506 

Material 

2024-T4 Aluminum 

2024-T4 

Elastic Modulus E (psi) 

io 7 

io 7 

PP Yield Stress a (psi) 
o 

50,000 

50,000 

Fragment Data 



Type 

T-58 Single Blade 

T-58 Single Blade 

Material 

SEL-15 

SEL-15 

Outer Radius (in) 

7.0 

7.0 

Fragment Centroid from Center of 
Rotation (in) 

4.812 

4.812 

Fragment Tip Clearance from Ring (in) 

1.657 

1.658 

Fragment Length (in) 

3.5 

3.5 

Fragment Length from CG to Tip (in) 

2.188 

2.188 

Fragment Weight (lbs) 

0.084 

0.084 

Fragment Moment of Inertia about its 
CG (in lb sec 2 ) 

2.163X10 -4 

2. 163xl0~ 4 

Failure Speed (RPM) 

15,374.3 

15,644.4 

Fragment Tip Velocity (ips) 

11,270. 

11,467. 

Fragment Centroidal Velocity (ips) 

7,748. 

7,884. 

Fragment Initial Angular Velocity (rad/sec) 1,610. 

1,638.3 

Fragment Translation KE (in lb) 

6,525. 

6,756. 

Fragment Rotational KE (in lb) 

280.4 

290.3 
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SUMMARY OP CIVM CALCULATION CASES FOR SINGLE-BLADE FRAGMENT IMPACT 
AGAINST A FREE COMPLETE CIRCULAR RING 
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SUMMARY OF CIVM CALCULATION CASES FOR SINGLF-RLADE FRAGMENT IMPACT 
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TABLE 4 


BLADE MOTION AND IMPACT DATA FOR CASE CR-5B(EL-PP) 
VS. CASE CR-7B (EL-PP-SR) 



Case CR-5B 




Case 

CR-7B 




(EL-PP) 




(EL-PP-SR) 


TAII 

ANI 

V 

V 

0)- 

ANI 

V 

V 

U). 

(Msec) 


y 

(in/sec) 

Z £ 

(in/sec) (rad/»ec) 


y 

(in/sec) 

z 

(in/sec) 

f 

(rad/sec) 

Pre-Impact 

0 

-7748. 

0 

+1610.0 

0 

-7748.0 

0 

+1610.0 

0 

1 

-6786.9 

-1125.3 

+2249.9 

1 

-6786.9 

-1125.3 

+2249.9 

16 

2 

-6604.2 

-1267.7 

+2266.5 

2 

-6523.8 

-1330.4 

+2273.9 

23 

5 

-6332.4 

-1476.1 

+2289.4 

6 

-6231.8 

-1551.7 

+2297.3 

27 

7 

-6216.2 

-1561.6 

+2297.3 

7 

-6099.5 

-1647.6 

+2305.8 

35 

10 

-6018.5 

-1702.0 

+2308.1 

11 

-5893.0 

-1792.2 

+2316.9 

36 

11 

-5871.2 

-1802.6 

+2314.3 

12 

-5877.9 

-1802.6 

+2317.6 

44 

— 

— 

— 

— 

15 

-5599.2 

-1989.5 

+2330.3 

45 

14 

-5614.4 

-1970.0 

+2322.0 

— 

— 

— 

— 

50 

15 

-5457.2 

-2068.6 

+2325.7 

17 

-5452.6 

-2084.6 

+2336.9 

61 

17 

-5300.2 

-2164.0 

+2329.4 

21 

-5132.4 

-2285.5 

+2354.6 

71 

19 

-4934.2 

-2377.3 

+2343.9 

__ 

— 

— 

-- 

72 

-- 

__ 

— 

— 

24 

-4933.6 

-2406.7 

+2373.7 

79 

— 

— 

— 

— 

25 

-4894.9 

-2430.2 

+2380.0 

85 

21 

-4781.7 

-2463.1 

+2360.5 

— 

— 

— 

— 

141 

— 

— 

— 

— 

26 

-3967.9 

-2765.6 

+2406.8 

144 

22 

-4449.1 

-2581.8 

+2370.7 

— 

— 

— 

— 

168 

25 

-3827.3 

-2786.4 

+2391.7 

— 

— 

— 

— 

178 

28 

-3724.7 

-2815.7 

+2396.5 

— 

— 

— 

— 

188 

— 

— 

— 

— 

28 

-3477.0 

-2895.3 

+2447.0 

189 

32 

-3569.7 

-2857.8 

+2405.6 

— 

— 

— 

— 

198 

35 

-3459.5 

-2886.0 

+2414.6 

— 

— 

— 

— 

210 

42 

-3353.0 

-2911.6 

+2426.9 

— 

-- 

— 

~ 

215 

— 

— 

— 


29 

-3453.0 

-2901.0 

+2451.3 

222 

45 

-3305.6 

-2922.4 

+2434.5 

— 

— 

— 

— 

246 

53 

-3233.8 

-2937.9 

+2450.2 

31 

-3287.9 

-2939.5 

+2488.1 

255 

58 

-3213.0 

-2942.2 

+2456.3 

— 

— 

— 

— 
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TABLE 4 (CONCLUDED) 



Case CR-5B 




Case CR 

-7B 




(EL-PP) 




(EL-PP- 

SR) 




V 

V 

0)^ 

ANI 

V 

V 


TAXI 

ANI 

y 

z 

f 


y 

z 

f 

(lisec) 


(in/sec) 

(In/sec) 

(rad/sec) 

ANI 

(in/sec) 

(in/sec) 

(rad/sec) 

264 

61 

-3199.8 

-2944.8 

+2460.8 

— 

— 

— 

— 

273 

63 

-3032.1 

-2971.7 

+2513.3 

— 

— 

— 

— 

274 

— 

— 

— 

— 

32 

-3110.8 

-2966.9 

+2544.9 

298 

65 

-2951.6 

-2982.0 

+2543.2 

39 

-2996.9 

-2982.2 

+2587.6 

307 

71 

-2888.0 

-2989.6 

+2567.9 

— 

— 

— 

— 

310 

— 

— 

— 

— 

42 

-2941.6 

-2988.9 

+2611.1 

316 

76 

-2841.8 

-2994.6 

+2587.4 

43 

-2939.2 

-2989.2 

+2612.2 

326 

80 

-2816.8 

-2997.1 

+2598.7 

46 

-2922.5 

-2991.0 

+2620.2 

327 

81 

-2812.7 

-2997.5 

+2600.6 

47 

-2912.0 

-2991.3 

+2621.5 

335 

84 

-2802.0 

-2998.5 

+2605.8 

II 

II 

II 

II 

345 

87 

-2794.9 

-2999.1 

+2609.5 

II 

II 

It 

II 

347 

88 

-2794.5 

-2999.1 

+2609.8 

ll 

II 

M 

II 

348 

89 

-2791.0 

-2999.4 

+2611.7 

II 

II 

•1 

It 

• 

tl 

It 

II 

M 

II 

II 

II 

It 


It 

II 

II 

II 

II 

II 

II 

II 


II 

II 

II 

II 

II 

II 

II 

II 

818 

89 

-2791.0 

-2999.4 

+2611.7 

47 

-2912.0 

-2991.3 

+2621.5 
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(a) Improved Method 


FIG. 4 FLOW CHART FOR SOLUTION PROCESS OF STRUCTURAL LARGE 
DEFLECTION ELASTIC-PLASTIC TRANSIENT RESPONSES 
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FIG. 5 NOMENCLATURE FOR GEOMETRY , COORDINATES , AND DISPLACEMENTS OF 
A CURVED-BEAM FINITE ELEMENT 
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FIG. 6 COMPARISON OF FINITE-ELEMENT WITH EXACT NORMAL MODE SOLUTION FOR THE SMALL-DEFLECTION 
RESPONSE OF A MECHANICALLY LOADED LINEAR-ELASTIC SIMPLY-SUPPORTED BEAM 




s 


where 

» effective shear carrying area of beam cross section 

w » transverse deflection 
M = bending moment 
S «« transverse shear force 

FIG. 7 GEOMETRY AND NONDIMENS IONAL QUANTITIES FOR A FREE-FREE 
BEAM SUBJECTED TO AN APPLIED CONCENTRATED TRANSIENT 
LOAD AT ITS MIDSPAN 
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(b) Clamped Ring 
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(b) Bounds of u) 

ItlctX 

FIG. 12 


2. 



MAXIMUM NATURAL FREQUENCY ( W xlO ) , RADIANS/SEC 

max , 



(a) The Maximum Natural Frequency 


FIG. 12 MAXIMUM NATURAL FREQUENCIES FOR SMALL VIBRATIONS OF 
THE CLAMPED BEAM 
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10 ELEMENTS PER HALF SPAN 
ELASTIC, PERFECTLY-PLASTIC BEAM 



(•Ni) Non.CH'i.jaa NvasaiN 
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TIME (MICROSECONDS) 

(a) Central-Difference Method 

FIG. 13 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED BEAM CALCULATED BY USING VARIOUS 
TEMPORAL FINITE-DIFFERENCE OPERATORS: CENTRAL-DIFFERENCE, HOUBOLT'S, AND NEWMARK'S 


10 ELEMENTS PER HALF SPAN 
ELASTIC, PERFECTLY-PLASTIC BEAM 
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(b) Houbolt’s Method 


10 ELEMENTS PER HALF SPAN 
ELASTIC, PERFECTLY-PLASTIC BEAM 
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TIME (MICROSECONDS) 







10 ELEMENTS PER HALF SPAN 
ELASTIC ,PERFECTLY-PLASTIC BEAM 







10 ELEMENTS PER HALF SPAN 
ELASTIC, PERFECTLY-PLASTIC BEAM 
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FIG. 15 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED BEAM CALCULATED BY THE 
CONVENTIONAL AND THE IMPROVED FORMULATION OF DYNAMIC EQUILIBRIUM 
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FIG. 16 COMPARISON OF DYNAMIC RESPONSES OF THE FREE RING OBTAINED BY USING 
CC VERSUS LC ASSUMED DISPLACEMENT FUNCTIONS 
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TIME (MICROSECONDS) 



FIG. 16 CONCLUDED 
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FIG. 17 COMPARISON OF DYNAMIC RESPONSES OF THE FREE RING (INCLUDING STRAIN-RATE EFFECT) 
OBTAINED BY USING LC VERSUS CC ASSUMED DISPLACEMENT FUNCTIONS 



(b) Outer Surface Strain Responses 
FIG. 17. CONCLUDED 
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FIG. 18 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED RING CALCULATED BY USING 
CC VERSUS LC ASSUMED DISPLACEMENT FUNCTIONS 



FREQUENCY (W x 10 ) RADIANS/SEC 



(a) MAXIMUM NATURAL FREQUENCY 

FIG. 19 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED BEAM CALCULATED 
BY USING THE LUMPED MASS MATRIX AND THE CONSISTENT MASS MATRIX 
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FIG. 20 CONTINUED 
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LUMPED MASS MATRIX 
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PIG. 23 EFFECTS OF INITIAL NODAL VELOCITIES AND DISCRETE ELEMENT SIZE ON DYNAMIC 
RESPONSES OF THE FREE RING 
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TIME (MICROSECONDS) 

FIG. 24 EFFECTS OF NUMBER OF SPANWISE GAUSSIAN INTEGRATION STATIONS ON DYNAMIC 
RESPONSES OF TI1E CLAMPED BEAM 




PIG. 25 EFFECTS OF NUMBER OF DEPTHWISE GAUSSIAN INTEGRATION POINTS ON DYNAMIC 
RESPONSES OF THE CLAMPED BEAM 
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THEORY : 1 9ELEMENTS/HALF 
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FIG. 28 EFFECTS OF STRAIN-RATE ON DYNAMIC RESPONSES OF THE CLAMPED RING 


MAXIMUM NATURAL FREQUENCY <U X 10 ) , RADIANS/SEC 



1 2 4 6 8 10 12 14 16 

NUMBER OF ELEMENTS PER HALF SPAN 

FIG. 29 COMPARISON OF THE MAXIMUM NATURAL FREQUENCY OF THE CLAMPED BEAM 
CALCULATED BY USING TIMOSHENKO-TYPE ELEMENTS VERSUS BERNOULLI - 
EULER-TYPE ELEMENTS 
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MAXIMUM NATURAL FREQUENCY (U) x 10 ) , RADIANS/SEC 

max 


3.0 


2.0 


2 3 6 9 12 15 18 

NUMBER OF ELEMENTS PER HALF RING 



FIG. 30 COMPARISON OF THE MAXIMUM NATURAL FREQUENCY OF THE FREE 
RING CALCULATED BY USING TIMOSHENKO-TYPE ELEMENTS VERSUS 
BERNOULLI -EULER-TYPE ELEMENTS 


244 




245 


500 1000 

TIMEC MICROSECONDS.) 

FIG. 31 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED BEAM CALCULATED BY USING THE 
TIMOSHENKO-TYPE ELEMENT VERSUS THE BERNOULLI-EULER-TYPE ELEMENT 
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FIG. 32 COMPARISON OF DYNAMIC RESPONSES OF THE CLAMPED RING CALCULATED BY USING 
TIMOSHENKO-TYPE ELEMENTS VERSUS BERNOULLI -EULER-TYPE ELEMENTS 




(*NI) NOIiVHYdHS aNVIdCIIW 3NriH3XN30 


247 


EXPERIMENT 
THEORY: EL-SH-SR 




(b) Deformation Profiles 


FIG. 33 CONTINUED 
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FIG. 33 CONTINUED 
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Surface Strain Histories 


EXPERIMENT 
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FIG. 33 CONCLUDED 
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FIG. 34a CONVERGENCE BEHAVIOR OF FE-PREDICTED PEAK RELATIVE CENTERLINE 
DEFORMATION FOR THE FREE-RING EXAMPLE 
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FIG. 34b APPROXIMATE CONVERGENCE BEHAVIOR OF FE-PREDICTED PEAK 
CIRCUMFERENTIAL STRAINS FOR THE FREE-RING EXAMPLE 



FREE RING 

FINITE-DIFFERENCE PREDICTIONS 
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FTG. 35b APPROXIMATE CONVERGENCE BEHAVIOR OF FD-PREDICTED PEAK 
CIRCUMFERENTIAL STRAINS FOR THE FREE-RING EXAMPLE 
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ROTATION 




MULTI -BLADED 
ROTOR 


(a) Single-Blade Fragment 




DISK 


(b) Multi-Bladed Disk Fragment 



(c) Multiple-Blade Fragments 

FIG. 37 SCHEMATICS OF SOME TYPES OF FAILED-ROTOR FRAGMENTS 
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NOTE: 1. U_ represents only the component of fragraent-c.g. velocity which 
is perpendicular to the (idealized-as-straight) ring segment 
imme diately prior to impact. 

2. U and U are, similarly, the perpendicular-direction components 
* of the ring segment nodal velocities. 

FIG. 39 SCHEMATICS OF FR&GMENT-RING COLLISION MODELS 
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NODE 2 


NODE 1 



(a) Consistent-Mass Model 



(b) Lumped -Mass Model 


FIG. 40 EXPLODED SCHEMATICS OF IDEALIZED COLLISION MODELS 
AT THE INSTANT OF IMPACT 







(a) Free 


(b) Smooth Hinge Support 



(c) Hinge-Restrained witn Elastic Foundation 
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(d) Point Elastically Restrained plus Elastic Foundation 



FIG. 42 SCHEMATICS FOR FREE AND RESTRAINED PARTIAL RINGS 



(a) Projection Inspection 

FIG. 43 INSPECTION FOR DETERMINING A COLLISION OF THE FRAGMENT WITH THE RING 







FRAGMENT COLLIDES WITH RING 



FIG. 44 INFORMATION FLOW SCHEMATIC FOR PREDICTING RING 
AND FRAGMENT MOTIONS IN THE CIVM APPROACH 
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(a) Complete Ring (Containment Device) 



bl 


b2 


b3 


(b) Ring Quadrant (Deflection Device) 


FIG. 45 SCHEMATICS OF EXAMPLE PROBLEMS ANALYZED BY THE CIVM APPROACH 



CONSISTENT MASS COLLISION MODEL 
° 9 ELEMENT (RQ-1B) 


o IQ ELEMENT (RQ-2B) 

a 15 ELEMENT (RQ-3B) 



(d) t « 



FIG. 46 CONCLUDED 




LUMPED MASS COLLISION MODEL 

o 9 ELEMENT 

o 10 element 

A 15 ELEMENT 




FIG. 47 CONCLUDED 


PREDICTIONS FOR a » 1 AND AN EL-PP RING 

a CM WITH 5-5-5-S DISCRETIZATION (CR-1B) 

a CM WITH 9-6-6-6 DISCRETIZATION (CR-3B) 



(a) t = 570 psec 

FIG. 48 COMPARISON OF CM WITH LM COLLISION MODEL PREDICTIONS FOR 
THE COMPLETE FREE RING SUBJECTED TO SINGLE-BLADE IMPACT 
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PREDICTIONS FOR e ■ 1 AND AN EL-PP RING 


° CM WITH 5- 5- 5- 5 DISCRETIZATION (CR-lB) 

A 04 WIXH 9-6-6-6 DISCRETIZATION (CR-3B) 

o LM WITH 10-6-6-6 DISCRETIZATION (CR-5B) 



FIG. 48 CONCLUDED 


272 


LOCATION 



1 i I i L. i — — _i i I i L_ 

rO CVl — O T 

(iNaoaad) 3 4 nitois aviiNaaaawnoaio 


273 


FIG. 49 EFFECT OF THE COEFFICIENT OF RESTITUTION ON THE PREDICTED INNER-SURFACE AND 
OUTER-SURFACE CIRCUMFERENTIAL STRAINS OF A BLADE-IMPACTED FREE COMPLETE RING 


ACCUMULATED NUMBER OF IMPACTS, ANI 



FIG. 50 EFFECT OF THE COEFFICIENT OF RESTITUTION ON THE ACCUMULATED 
NUMBER OF IMPACTS OF A BLADE-IMPACTED FREE COMPLETE RING 
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EXPERIMENT 


X CASE CR-5B (EL-PP, e - 1} 

A CASE CR-7B (EL-PP-SR, e - 1) 

RING BEFORE INITIAL IMPACT 



(a) Time after Initial Impact (TAII) ■ ISO usee 

FIG. 51 COMPARISON OF PREDICTIONS WITH EXPERIMENT FOR THE FREE COMPLETE 
RING SUBJECTED TO SINGLE-BLADE IMPACT IN NAPTC TEST 88 
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o 


EXPERIMENT 


X CASE CR-5B (EL-PP , e - 1) 

a CASE CR-7B (EL-PP-SR, e - 1) 

RING BEFORE INITIAL IMPACT 



(b) TAII - 570 psec 
FIG. 51 CONTINUED 
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o 


EXPERIMENT 


x CASE CR-5B (EL-PP , e * I) 

A CASE CR-7B (EL-PP-SR, e « 1) 

RING BEFORE INITIAL IMPACT 



(c) TAII * 810 ysec 
FIG. 51 CONCLUDED 
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O EXPERIMENT 

X CASE CR-11B (EL-PP-SR, e * 0) 

A CASE CR-I0B (EL-PP-SR, e - 1) 

RING BEFORE INITIAL IMPACT 



FIG. 52 CONTINUED 
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o EXPERIMENT 

X CASE CR-I1B (EL-PP-SR, e - 0) 

a CASE CR-10B (EL-PP-SR, e - 1) 

RING BEFORE INITIAL IMPACT 



FIG. 52 CONCLUDED 
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EXPERIMENT (WITH REVISED INITIAL IMPACT INSTANT) 
CASE CR-10B (EL-PP-SR, e * 1) 

RING BEFORE INITIAL IMPACT 


0 



0 

(a) TAII * 570 usee 

FIG. 53 COMPARISON OF PREDICTIONS WITH EXPERIMENT FOR THE FREE COMPLETE 
RING SUBJECTED TO SINGLE-BLADE IMPACT IN NAPTC TEST 91, WITH A 
REVISED INSTANT OF INITIAL IMPACT 
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EXPERIMENT (WITH REVISED INITIAL IMPACT INSTANT) 


a — 


— CASE CR-10B (EL-PP-SR, e - 1) 

— RING BEFORE INITIAL IMPACT 



FIG. 53 CONCLUDED 
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(d) t 


IC END CASE HQ-5B 
HINGED END CASE RQ-8B 



FIG. 54 CONCLUDED 






APPENDIX A 


DESCRIPTION OF THE MECHANICAL SUBLAYER MODEL FOR STRAIN-HARDENING, 

STRAIN-RATE SENSITIVE MATERIAL BEHAVIOR 

As discussed in Subsection 2.3.2, the yield surface of certain materi- 
als will change in case of continued straining beyond the initial yield. The 
change of the yield surface that characterizes the strain-hardening behavior 
of the material depends on the loading history. In the present analysis, the 
strain-hardening behavior of the material is accounted for by using the 
"mechanical sublayer model" (Ref. 20, 21, 46, and 48). In order that the 
present report be reasonably self-contained, the mechanical sublayer model is 
described in this appendix. 

In the mechanical sublayer model, the uniaxial tension (or compression) 
test stress-strain curve of the material is first approximated by (n+1) piece- 
wise-linear segments which are defined at coordinates [ (O^, £^) , k = 1, 2, ...n], 
as depicted in Fig. A. la. Next, the material is envisioned as consisting, at any 
point in the material, of n equally-strained "sublayers" of elastic, perfectly- 
plastic material, with each sublayer having the same elastic modulus E, but an 
appropriately different yield stress (see Fig. A. lb) . For example, the yield 
stress of the kth sublayer is 

(Tok = E , ( k = i, 2, , n ) (A.i) 

Then, the stress value, associated with the kth sublayer can be defined 
uniquely by the strain history and the value of strain and strain-rate presu.it 
at that point. Taken collectively with an appropriate weighting factor for 
each sublayer , the stress , 0 , at that point corresponding to strain £ may be 
expressed as 

O' ~ Ck Kit) (A. 2) 

k = 1 

where the weighting factor for the kth sublayer may readily be confirmed to 
be 

v -'fc — fr (a. 3) 


a - Ek-i 
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The elastic perfectly-plastic and linear strain-hardening constitutive rela- 
tion may be treated as special cases. In the case of elastic perfectly- 
plastic behavior, there is only one sublayer, and in the case of linear strain- 
hardening material there are two sublayers and the yield limit of the second 
sublayer is taken sufficiently high so that the deformation in that sublayer 
remains elastic. 


From the computational point of view, the use of the mechanical sub- 
layer model is very convenient to analyze problems with general loading paths 
including loading, unloading, reloading, and cyclic loading. Its features 
include the "kinematic hardening rule" which takes the Baushinger effect into 
account (see Fig. A.lc). Also, this mechanical sublayer model may readily 
accommodate the strain-rate effect. Figure A. 2a illustrates schematically the 
uniaxial stress-strain behavior for a strain-rate dependent, elastic, perfectly- 
plastic material whose rate dependence is described by Eq. 2.87, 

= S ( 1 + ' 

while Fig. A. 2b depicts the corresponding behavior for a strain-hardening ma- 
terial which is represented by the mechanical sublayer model, each sublayer of 
which has the same values for the strain-rate constants D and p. For this 
special type of rate-dependent strain-harde.iing material, the stress-strain 
curve at a given strain rate e is simply a constant magnification of the static 
stress-strain curves along rays emanating from the origin. Ho - . ’ever , for strain- 
hardening material whose strain-rate behavior is not one of simple magnifica- 
tion, the strain- rate behavior can often be approximated adequately by employ- 
ing appropriately different values of D and p for each sublayer; the resulting 
behavior is shown schematically in Fig. A. 2c. 

One may generalize this uniaxial behavior to the two- or three- 
dimensional stress case by adopting, for example , the Mises-Hencky yield con- 
dition, Eq. 2.79, and flow rule, Eq. 2.82, and applying them to each sublayer 


t ) 

/ (2.87) 
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of the mechanical sublayer model. The strain-rate dependence may also be gen- 
eralized by assuming that the e of the one-dimensional case may be replaced 
by the second invariant of the deviatoric strain-rate tensor as defined by 


3 r 7 * 




(3.64) 


In terms of the finite increments Ay^ of strain determined in each 
timewise calculation step of the present procedure, the "replacement £" 
given by Eg. 3.64 becomes: 

i 

| ( ^ -t , vJ I , .-'Y .2 


£ = 


-at 


It Ar ’-T (Ar £ f ] 


(A. 4) 
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UNIAXIAL STRESS, a UNIAXIAL STRESS, 


(<r 4 <4) 



STRAIN , e 

(b) PROPERTIES OF THE ELASTIC, PERFECTLY- PuA'VTiC SUBLAYERS 


FIG. A.l APPROXIMATION OF A UNIAXIAL STRESS-STRAIN CURVE 
BY THE MECHANICAL SUBLAYER MODEL 
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STRAIN, € 

(b) SPECIAL STRAIN HARDENING MATERIAL 

FIG. A. 2 SCHEMATIC OF STRAIN-RATE DEPENDENT UNIAXIAL 
STRESS-STRAIN CURVES 
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APPENDIX B 


DEFINITION OF FINITE ELEMENT MATRICES INDICATED 
SYMBOLICALLY IN THE TEXT 


The various matrices which are used symbolically in Section 4 for an 
arbitrarily curved beam element are documented in this appendix, and their 
specialization to represent simple circular ring and straight beam elements 
are also listed. 


B.l Bernoulli-Euler Type Deformation Behavior 

B.1.1 Bernoulli-Euler Type Curved Beam Element with 
CC Assumed-Displacement Functions* 


iij - ' 

w J 

-(U'f'j 
2 *8 ' 

II Pi 

8*1 




(4. 
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' — r~ 

II 
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Vi 

w, V; V 
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'?<• 

vV t > 1 ^>1 
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0 w 

(4. 

.17) 
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1*8 

fj) + l 


{d,J 

5 « t 1*8 

h\ 


(4. 

.24) 

Kr — L 

C>3J 

1 < 8 

1 if) 





(4. 

.24) 

(m) = j 

f 1i~ 1 

( 

B] { 

N(^) ^ 



(4. 

.31b) 

J 

u 

8 <-3 

3 * 3 

3 



If) = 

[?, 
1 1< 

’ [ N <l>fj 

r£i 

Hi 




(4, 

.34a) 

II 

fli* 1 

4< 

, | a j L ■* 

l°3 

j M ) ^ ? 



(4, 

.37a) 

f h) = 

rU 

t a) l d, 

- J L 

dy l 



(4, 

.37a) 


* N ^ (4.40a) 

f 

Equations not numbered with a prefix B refer to pertinent equations in the text. 
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{^1 = "</*7 EU( r i f j fAMj if} > f 

,h- i (4.40d; 

+lf£kh I fy + r ifj f A )liu ff f )[I>4 i4j '? If i> 

I ff j = f/ /J‘"( E e f j A } + 5 E e*l d,\ ) i V t 

I f Pi = (} (!" E^IdAlDJ d V, 4) ( 4 . 40cl 

A ic ’ 

In the following, the matrices fU^], fA] ,UQj, x = ^ 2 . 3) , (Nf^J.and 
[Bj appearing in the above expressions are listed for an arbitrarily-curved 
beam element (see Fig. 5): 


C ° S ^ ~Zcos<p+ysiN<P % 0 0 *1* >[ 3 

—sin <P cos 4* z sim 4>i-y cost o y 2 y 3 o o (4ml3) 


COS^ SI Mb Q 0 

-SIM4>£ <L0S<Pi o 0 

o 0 I 0 

, SI Mb -costz 

M) = R* R* 0 1 

cos4h jtw4, X +| S'N^ ( {?*, 

-vJINi, C0S4^ 0 

0 0 I 

SI N 4^1 -COS& I _y C0s4>+ , 

/"S£. Riv. R,*, I 

lD,j=lO 00 I n y R f/ R 
L ^j-LO O / - ’/ft 2? 3? 2 

lC,j=i° o o A-?lf -<? 


0 0 


o o 


0 0 %, t 

? ( vi 7,., o o 

X, X 

o o 27 ,,, 37 ,'. 

2 ? 3 7 ' J f jA"] ( 8 . 2 ) 

-l‘/R -f/Rj[A"] a . 3> 

K‘] (B .4, 
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[N<?']= 


[B] 


COS# N 4" 

cos 4 s 

0 0 


"l _ 0 
m 

0 7 


-2 cos4+ysin4‘ ^ 0 0 Y 


zsiNQ+ycosf o *1“ Y 

I -k 2 ? 3 ? 


•A 

R 


V 

0 

-t 

R 


(A"J 


(B. 5) 


(B.6) 


where 


m = mass per unit length of the beam element 
I = mass moment of inertia of the beam cross section 


Note that, in the above equations, only the matrix {A ] is given. In 
practice, one forms its inverse ,[A '] , by means of available standard computer 
subroutines. Also note that subscript i ( or i+1 ) is used to denote the value 
of the quantities at node i ( or node i+1 ) . 

For application to simple circular ring elements, the geometry and 
nomenclature of a typical circular ring element is shown in Fig. B.l where the 
local and global coordinates are arranged to take advantage of the symmetry of 
the ring element's geometry. The matrices ( U , l A) , (. L A.' 1 , < — I, ^ , 3) 

[f\j ( ^)] [m] /and f f<J for a simple circular ring element are listed in the 
following : 
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(B.8) 
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(B. 9 ) 
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where, in Eqs. B.13 and B.14 , I = $ kh. k, = EI;k^ and k 2 ' 

f 0 = mass density per unit volume of the unde formed beam 
{7 =the width of the beam 
k =the thickness of the beam 
E "Young’s modulus 


(B.14) 


eU 3 

12 
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B.1.2 I -AouXli-Euler Type Curved Beam Element with 
LC Assumed-Displacement Functions 


/• 




tfe1=/£ I 


(B.X5) 

(B.16) 

(B„17) 

(B.18) 

(B.19) 

(B . 20) 
(B.21) 

(B.22) 

(B.23) 


iff j = (( {}"( EE^I^j + 5 ^3 j ) d V'n (B.25) 

i ff \ = f( fj' E E.* f L D zJ d V„ f j\ (b.26) 

r ’a u 


In the folowing , the matrices [(_! C^j] ~ -1,2, 3 ) and 

f B ] appearing in the above expressions are listed for an arbitrarily-curved 
beam element ( see Fig. 5 ) : 
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For application to a circular ring element ( see Fig. B.l ) , the matrices 

IA 3, (lAiJ,*'=V. 3), [N^O, [m] ,and[kl sure listed in the 

following: 
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The corresponding quantities for a straight beam element (see Fig. S.2 ) 
can be obtained by setting -^ = O, A = 0 , and A = d'j in the preceeding 
equations. Thus, one obtains 
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B.2 Timoshenko-Type Curved Beam Element 

For the Timoshenko-type curved beam element, only the element properties 
with the linear displacement interpolation functions are presented here. The 
element properties with higher order assumed displacement functions may be 
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derived in a similar manner as described in Subsection 4.4 . 
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In the following, the matrices [U^J, £^1, ( L^J, i- /, 2 , 3,4) fM ( ? J 3, 
and [B] appearing in the above expressions are listed for an arbitrarily- 


curved beam element 
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For application to a circular ring element ( see Fig. B.l ) , the matrices 
[U<7 ( l4j, i. 2,3,4), [m] ' and rk] sure listed in the 

following: 
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FIG* B.l SIMPLE CIRCULAR RING ELEMENT 


z,w 



FIG. B.2 STRAIGHT BEAM ELEMENT 
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APPENDIX C 

MIXED FINITE-ELEMENT MODEL BASED UPON REISSNER'S VARIATIONAL PRINCIPLE 


This appendix contains, for the purpose of demonstrating the use of 
different variational principles, the formulation of the equation of motion 
based on a mixed finite-element model: (1) by assuming a displacement field 
which is continuous over the entire solid, (2) by assuming a stress field for 
each individual element. However, no evaluation of this model is made. 

In the Principle of Virtual Work, the independent quantity which is 
subject to variation is the displacement; in Reissner's variational principle, 
the independent quantities which are subject to variation are stress and 
displacement. Reissner's variational principle may be stated mathematically 
as follows (Ref. 159} : 

Sir R ={({S(s‘ i y, i )*v.-lB-Sw*Si-([ k-£) St 'aa. =o 

v. A" K- 1 ' 

where 

SB =((( x f SS‘M=f(f ( Cjjia S k + Xf ) J 5‘W. 

= the variation of the complementary 
work done in an arbitrary infinitesi- 
mal virtual stress OS 1 - 1 increment 

/ (C.2) 


(C.3) 


~ l ( Vi >i + V * v % ) 

/\ov - portion of A q over which the displacement v^ is prescribed 
C,J kX = elastic compliance tensor 
and^VV d I ~f* Y.tare defined as previously stated in Section 2. 

I 


In applying Eq. C.l to the finite-element analysis for a solid continu- 
um which is conceptually subdivided into N discrete elements, one can write 
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^l{(f[S(S%)dV^SB n ^W n ^I n -ff< * - a *TdA,-C . ) 

V„ As, (c.4) 

where the term C arises from a possible jump function of the derivative of 
n 

across the interelement boundaries. If the assumed displacement functions 

are such that v. and v. . cure continuous across the interelement boundaries, 

x X,] 

then C = 0 (Ref. 7) . 
n 

The assumed approximate stress interpolation functions which need not 
satisfy the equilibrium condition for each individual element are expressed 
using a finite number of independent parameters {a}: 

S* ( = L R‘ ’ J K (C.5) 

and the displacements are approximated by interpolation functions [N] to 
satisfy the interelement displacement and slope compatibility conditions 
anchored to nodal generalized coordinates {q}: 

l v j = [N]ffJ 

By applying Eq. C.6, Eq. C.3 becomes 

Y*y = L -1 \ f\ + I L f J L D ?' J \f\ (c . 7) 

Using Eqs. C.5, C.G, and C.7, and if the boundary displacement can be 
made to satisfy the prescribed value, Eq. C.4 may be expressed in the form 
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where the following quantities are evaluated for each finite element: 

ir}=(([ f R ,? ] C,,u lR“j 

v. 


(C.8a) 
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[*} = /f/ i^ 7 j ( i- % J ^ J f D *<\ t D f J 'fr\)JV n (c - 8c) 

fAf} ( ^ (C. 8d) 


312 



[f > ]=[({ I j i D )'J (lR"j !<*\UV„ 


ff| =fff (Nfp.lfWK, *(( CNafflHA 


(C.8f ) 


r^j = ((( [n]>. [N]dv„ 


(C.8g) 


The quantities {a} can be varied independently for each finite element; 
hence, from Eq. C.8 it may be concluded that 


r r) iofj - f s j - jt) = o 


Mj = t K ( i sj + !*) ) 


and, hence, it follows from Eqs. C.9 and C.8 that 


(C.9a) 


X Lip ( lfl + [hj{ff}-t-fj- + [M}tjf.j)=0 {c . 10) 


It is seen that Eq. C.10 is of the same form as Eq. 3.14 which is 
associated with the assumed-displacement type of finite-element model. Thus, 
transforming the element nodal displacements {q} to independent global dis- 
placements {q*} of the discrete-element assembly as described previously will 
lead to the same form of matrix expression as given by Eq. 3.17 which represents 
the "improved formulation" form of the dynamic equilibrium equations: 

rvuffc'j + {p} + cm ffj = r f*i 

However, it should be noted that the numerical values of the terms in [M] , (p), 
[H] , and {F*} of Eq. C.ll can be, in general, different from the corresponding 
symbolic quantities in Eq. 3.17. 

Given a set of initial conditions {q*} , {q*} and {F*} at time zero, 

o o o 

and the proper boundary conditions, the system of second-order differential 
equations, Eq. C.ll, together with Eq. C.9 may be solved in a step-by-step 
£i.xn&wxs@ fashion. Lat it h© as sinned, that at a typicaX time instant one 
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knows the generalized nodal displacement the strain component , 

and the plastic strain component at any station in each finite 

element. The stress parameter {a} for each individual finite element may 

in 

be calculated from Eq. C.9, where {s} and {t} can be evaluated from Eqs. C.8b 

m m 

and C.8c, respectively. Having the stress parameter {a} at time t , one can 

m 

calculate {p} and [h] for each finite element by using Eqs. C.8d and C.8e, 
m m 

respectively. Now, assemble this information for the discretized structure. 

Then Eq. C.ll together with the use of a proper timewise finite-difference 

operator to approximate {q}^ can be employed to determine the displacements 

{q*} or the displacement increments {Aq*} at the next time instant t since 

all quantities except {q*} in Eq. C.ll are known, and hence {q*} can be readily 

m m 

calculated. Having the displacements or displacement increments at time t , , the 

m+i 

strain (or strain increment) at any station in each element may be found from Eq. 

C.7. With the displacement and strain available, the desired stress parameter at 

time instant t , can be determined through the use of proper elastic-plastic con- 
m+l 

stitutive relations, and Eqs. C.5 and C.9. Then Eq. C.ll can furnish the displace- 
ments for the next time step. The process is cyclic thereafter. 

It should be noted that the application of the mixed finite-element model in 
this form does not provide any particular advantage over the assumed-displace- 
ment finite-element model since the interelement boundary compatibility is 
still required. However, the use of the mixed finite-element model has its 
merit in the elastic analysis of plates and shells (Ref. 7) . But for the 
elastic-plastic analysis, the proper interaction (yielding) surface between 
stress and moment resultants has not been found; such information, however, would 
be needed in this mixed method. 

Finally, it can be expected that assumed-stress finite-element models, 
such as the assumed-stress hybrid model, the equilibrium model, etc., (Ref. 7) 
will not be as accurate as the assumed displacement finite-element model for 
elastic-plastic transient structural analyses, unless the time-step size is 
made sufficiently small, since the stress-strain curves for many structural 
materials are usually very flat in the plastic range; thus, a small error in 
the strain will produce a small error in the stress, but on the other hand, a 
small error in the stress will result in a much larger error in the strain. 
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